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Abstract 


A  series  of  lectures  are  presented  on  the  topic  of  the  location  and  identification  of 
compact  objects  by  low  frequency  electromagnetics.  These  lectures  were  presented 
as  a  portion  of  two  graduate  level  courses  in  electrical  engineering  at  the  University 
of  Toronto  in  1985  and  1987.  Magnetostatics,  electrostatics  and  electromagnetic 
induction  techniques  are  discussed  in  detail. 
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In  the  Fall  of  1985,  a  pilot  course  was  offered  to  graduate  students  in  electrical  engineering  at  the 
University  of  Toronto.  The  course,  entitled  “Topics  in  Electromagnetic  Theory  :  Latest  Techniques 
in  Remote  Sensing”,  was  taught  by  Professors  J.L.  Yen  and  K.  Iisuka  of  the  Department  of  Electrical 
Engineering  and  me.  Professor  Yen  lectured  on  microwave  radiometry,  Professor  Iisuka  covered  subsurface 
rada  and  I  dealt  with  low  frequency  electromagnetics. 

The  material  met  with  considerable  enthusiasm  to  the  end  that  it  was  offered  again  in  the  Fall  of  1987 
under  the  new  title  “Electromagnetic  Remote  Sensing"  with  the  same  lecturers  and  general  topics. 

These  notes  are  a  compilation  of  the  material  that  V7as  presented  for  the  low  frequency  electromagnetic 
remote  sensing  portion  of  both  courses.  The  general  subject  has  been  covered  extensively  in  textbooks  in 
the  context  of  geological  exploration.  Objects  of  interest  are  typically  considered  to  be  very  large,  often 
infinite  in  extent  in  one  or  two  dimensions.  I  have  approached  these  lectures  from  the  viewpoint  of  the 
location  and  identification  of  relatively  small  compact  objects  by  low  frequency  electromagnetics.  This  is 
an  area  which  has  been  surprisingly  neglected  in  textbooks  even  though  it  is  a  very  relevant  topic  which 
is  of  practical  interest  to  many  researchers  in  areas  as  diverse  as  submarine  detection,  pipe  location, 
archaeology  and  mine  and  artillery  shell  detection.  I  hope  these  notes  will  All  that  gap  in  the  available 
textbooks. 

Finally,  I  must  thank  my  colleague  Dr.  Yogadhish  Das.  Much  of  the  work  presented  in  these  notes  rep¬ 
resents  published  (and  some  unpublished)  collaborative  research  by  the  two  of  us.  Our  general  discussions 
have  been  invaluable  and  have  had  a  distinct  influence  on  these  notes. 


1  Low  Frequency  Electromagnetic 
Remote  Sensing 

The  present  course  is  concerned  with  remote  sensing,  particularly  as  it  relates  to  electromagnetics.  Remote 
sensing  may  be  defined  as  the  acquisition  of  information  about  an  object  without  physical  contact. 
Although  a  host  of  sensors  might  qualify  as  remote  sensors  under  this  broad  definition,  remote  sensing 
usually  refers  to  the  collection  and  processing  of  information  about  the  earth  or  objects  on  the  earth 
through  the  use  of  photographs  and  sensor  data  acquired  from  a  satellite  or  aircraft.  We  will  refer  to  this 
as  “classical  remote  sensing  (CRS)”. 

The  problem  of  detecting  objects  concealed  from  view  at  distances  of  between,  say,  1  and  100  times 
a  characteristic  dimension  of  the  object  is  an  area  of  remote  sensing  which  is  of  much  interest.  We  shall 
call  this  problem,  for  lack  of  a  better  term,  “quasi-  remote  sensing  (QRS)”.  The  boundary  between  QRS 
and  CRS  is  fussy.  Computed  axial  tomography  might  be  considered  to  QRS  when  used  to  image  the 
body  or  a  buried  object  in  soil,  but  tomography  using  radar  signals  to  image  an  airplane  lies  more  in 
the  realm  of  CRS.  Detecting  ore  bodies  using  an  airborne  electromagnetic  induction  detector  might  be 
considered  CRS  but  locating  the  same  bodies  from  the  ground  using  the  same  type  of  detector  might  be 
QRS.  We  shall  not  attempt  to  clearly  delineate  the  boundary  between  CRS  and  QRS  but  rather  shall 
give  representative  examples  of  the  latter.  Some  of  these  are: 

1.  Geophysical  exploration  of  compact  ore  bodies  from  the  ground  (magnetics,  electromagnetic  induc¬ 
tion,  electrostatics) 

2.  Nondestructive  testing  (electromagnetic  induction,  neutron  tomography,  x-ray  imaging) 

3.  Detection  of  submarines  from  aircraft  or  ships  (magnetics,  electromagnetic  induction) 

4.  Detection  of  land  vehicles  using  buried  sensors  (magnetics,  electromagnetic  induction) 

5.  Detection  of  buried  ordnance  (artillery  shells  and  mines)  (low  frequency  electromagnetics) 

6.  Detection  of  trapped  miners  (low  frequency  electromagnetics) 

7.  Imaging  of  body  parts  (x-rays,  electrostatics,  acoustics) 

8.  Measurement  of  electric  fields  of  body  organs  (electrocardiogram,  electroencephalogram) 

8.  Measurement  of  maguclic  fields  of  body  organs  (magnetocardiogram,  magnetoencephalogram) 

10.  Measurement  of  contaminants  in  the  body  such  as  dust  in  lungs  or  metal  fragments  (magnetics, 
electromagnetic  induction) 

In  the  previous  list,  some  technologies  applicable  to  each  problem  have  been  noted,  although  this  is 
by  no  means  an  exhaustive  compilation.  To  illustrate  the  wide  spectrum  of  methods  available  fot  QRS, 
we  shall  focus  for  a  moment  on  the  problem  of  delecting  unexploded  ordnance.  Table  1.1  shows  the  main 
characteristics  of  the  three  major  types  of  unexploded  ordnance.  Keeping  these  traits  in  mind,  one  can 
think  of  a  number  of  possible  methods  which  could  conceivably  detect  buried  ordnance.  These  arc  shown 
in  Table  1.2,  where  they  are  categorised  as  “Q"  or  “R"  for  quasi-remote  or  remote  sensing. 

A  discussion  of  the  principles  and  techniques  behind  all  the  methods  listed  is  beyond  the  scope  of  the 
course,  however,  details  may  be  found  in  (22). 

In  keeping  with  the  theme  of  the  course,  we  shall  restrict  ourselves  to  quasi-remote  sensing  methods 
which  ate  related  to  electromagnetism.  Although  one  could  argue  eloquently  that  all  the  above  methods 
are  indirectly  related  to  electromagnetism,  we  shall  be  less  philosophical  and  focus  on  those  that  directly 
involve  the  generation  and/or  measurement  of  electromagnetic  fields.  These  include  magnetostatics, 
electrostatics,  electromagnetic  induction,  ground  penetrating  radar,  magnetic  resonance,  and  the  sensing 
of  ultraviolet,  visible  and  thermal  infrared  radiation.  The  last  three  are  usually  considered  to  be  classical 
remote  sensing  techniques  and  ground  penetrating  radar  is  covered  elsewhere  in  the  course.  Although 
magnetic  resonance  involves  the  generation  and  tensing  of  radio  frequency  radiation  coupled  with  static 
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magnetic  fields,  it  requites  a  quantum  mechanical  treatment  and  could  easily  fill  one  or  more  courses  by 
itself. 

Thus,  we  will  restrict  this  portion  of  the  course  to  magnetostatics,  electrostatics  and  electromagnetic 
induction,  particularly  as  they  apply  to  detecting  nearby  concealed  objects.  In  particular,  we  shall 
concentrate  on  “compact”  objects,  that  is,  those  which  have  a  finite,  closed  bounding  surface  and  hence 
are  finite  in  extent.  All  three  methods  are  in  the  low  frequency  portion  of  the  EM  spectrum  (DC  to  sub- 
RF).  This  is  necessary  because  the  materials  concealing  the  object  of  interest  attenuate  electromagnetic 
radiation  and  the  attenuation  generally  increases  with  frequency.  For  example,  Table  1.3  shows  the  skin 
depth  at  a  few  frequencies  for  two  of  the  most  common  barriers  encountered  in  practice,  saline  (the  major 
EM  attenuating  constituent  of  the  body)  and  soil. 

Since  minimum  torso  thickness  is  approximately  0.16  m  and  minimum  soil  penetration  depths  are  of 
the  order  of  1  to  2  m,  solely  from  the  standpoint  of  signal  attenuation  it  is  desirable  to  use  frequencies 
less  than  1  MHs.  There  is,  however,  a  tradeoff.  Because  of  the  small  distances  and  long  wavelengths, 
measurements  are  usually  in  the  near  fieid.  Thus,  the  approximations  of  geometrical  or  physical  optics 
do  not  hold  and  imaging  of  a  source  is  not  a  straightforward  matter.  In  fact,  as  we  shall  see,  the  inverse 
problem  at  these  low  frequencies  is  quite  challenging. 

2  Magnetostatic  Methods 

2.1  Introduction 

The  object  of  magnetostatic  methods  is  to  measure  the  static  magnetic  field  associated  with  a  ferrous 
object  and  to  use  this  information  to  locate  and  identity  the  object.1  Ferrous  objects,  when  placed 
in  a  magnetic  field,  such  as  that  of  the  earth,  acquire  magnetisation  of  their  own.  This  is  normally 
referred  to  as  “induced  magnetisation”.  The  magnitude  and  direction  of  the  induced  magnetisation  is 
a  function  of  the  ferromagnetic  susceptibility  of  the  body,  its  shape  and  its  orientation  with  respect  to 
the  ambient  field.  Ferrous  objects  may  also  have  permanent  or  intrinsic  magnetisation,  usually  called 
“remnant  magnetisation".  Remnant  magnetisation  is  a  function  solely  of  the  thermal,  mechanical  and 
geomagnetic  history  of  the  object. 

We  will  be  dealing  primarily  with  "compact"  objects.  The  spatial  extent  of  the  measurcable  fields 
associated  with  such  objects  is  generally  less  than  a  few  meters  and  the  time  taken  to  measure  the  fields 
is  leu  than  a  minute.  U  is  important  that  the  earth's  field  not  change  significantly  during  that  period. 

2.2  Geomagnetic  Field 

When  the  geographical  variations  are  a  vet  aged  out,  the  earth's  magnetic  field  is  very  similar  to  that  of 
a  magnetic  dipole  located  ~  400  km  (*v  X  0f  the  earth's  radius)  due  north  of  the  geometrical  center  of 
the  earth  (Fig.  2.1,  Table  2.1).  Its  axis  lies  in  the  meridian  plane  defined  by  60°  W  and  it  lilted  11° 
from  the  geographical  N-S  axis.  The  south  pole  of  the  dipole  points  toward  geographical  north.  The 
magnitude  of  the  dipole  moment  is  approximately  8  x  10**  Am*  which  produces  an  average  polar  field 
of  65000  nanoTeslas  (nT)  and  on  average  equatorial  field  of  35000  nT  at  the  earth's  surface  (Figs.  2.2, 
2.3).  This  average  gcomsgnetic  field  is  mainly  due  to  currents  in  the  highly  conductive  core  of  the  earth 
(self-exciting  dynamo  action  of  thermal  currents) 

The  geomagnetic  field  possesses  a  high  degree  of  spatial  uniformity.*  For  example,  near  latitude 
45° N  where  F  is  approximately  45000  nT,  the  field  gradient  is  approximately  20  nT/km  in  altitude  and 
5  nT/km  in  latitude.8  If  the  sensor  is  leu  than  5  m  above  ground  surface,  iuhomogeneity  of  magnetic 
properties  of  soil  and  magnetotclluric  currents  may  disturb  the  local  '  ^ats-vl  uniformity.  Buried  ferrous 
objects,  geological  anomalies  and  pockets  of  residual  magnetisation  •'  V.i  la* Wr  sometimes  the  result 

1  Sialic  Magnetic  field*  are  mtuuitd  for  other  reason*  too.  For  tramp le  oat  ca*  tludy  Ike  Add ItxU u  la  interplanetary 
field  meMurenteati  or  deltnciae  ciimsl  diilribullctn*  a  it  Ida  Ike  hnmaa  body  from  biornagnttie  measurement*. 

*TkU  ie  because  large  variations  art  unlikely  over  Ike  large  volume*  of  fluid  la  Ike  earth'*  core  and  Ike  high  cole 
temperature  eradicate*  large  contrast*  afler  a  short  time. 

5  The  korUoalal  (radical  **  0  -  lOnT/lua  from  equator  lo  poiei  -0.03  aT/m  at  Magnetic  poieai  -0.015  bT/m  at  magnetic 
cqaatori  or  ±Q.047>’  aT/m  (+  north  of  equator,  *  tooth,  F  La  Oc). 
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of  currents  due  to  ancient  lightning  strikes)  will  also  disturb  the  uniformity  but  may  be  detectable  at 
standoffs  greater  than  5  m. 

Temporal  fluctuations  are  25%  attributable  to  currents  internal  to  the  earth  and  75%  due  to  external 
currents.  The  latter  are  mainly  located  in  the  ionosphere  and  variations  are  caused  primarily  by  fluctu¬ 
ations  in  the  solar  wind.  The  internal  and  external  currents  are  coupled  and  lead  to  a  complex  temporal 
variation  which  can  be  resolved  into  several  components  (Figs.  2.4-2.6): 

1.  Secular  variations:  Occur  at  a  rate  of  typically  30  nT/year. 

2.  Diurnal  variations:  At  mid-latitudes,  F  decreases  rapidly  soon  after  sunrise,  reaches  a  minimum 
at  noon  and  increases  during  afternoon  and  night.  The  amplitude  of  oscillation  is  a  maximum  in 
June  (~  ±25  nT)  and  is  minimum  in  January  (~  ±2.5  nT).  The  direction  of  the  field  also  changes 
by  a  similar  amount.4 

3.  Fast  variations  (micropulsations):  There  are  several  types  including  Eschenhagen  fluctuations 
(amplitude  ~  0.5  -  5.0  nT,  period  ~  25  sec,  total  duration  ~  few  minutes);  hydromagnetic  pearls 
(amplitude  ~  0.02  nT,  period  ~  0.3  -  3.  sec);  oscillation  of  the  earth/ionosphere  cavity  (amplitude 
~  0.002  nT,  period  ~  0.025  -  0.2  sec). 

4.  Exceptional  variations  (magnetic  storms):  These  happen  several  times  per  year  and  are 
linked  with  solar  activity.  In  fact,  they  occur  ~3Q  hr  after  strong  pulse-like  eruptions  on  the  solar 
surface  (solar  flares).  The  typical  variations  are  from  500  -  1500  nT  with  a  wide  range  of  frequency 
components.  In  extreme  cases,  the  fluctuations  are  so  strong  as  to  produce  visible  deflections  in 
compasses."  During  magnetic  storms,  the  only  thing  that  can  be  measured  are  magnetic  storms. 
Such  storms  may  last  several  hours. 

2.S  Interplanetary  Field 

Remote  tensing  also  includes  measurements  in  outer  space  and  magnetic  measurements  are  a  very  im¬ 
portant  component.  There  ore  three  main  areas  of  interest: 

1.  300km  <  r  <  800km  :  Measurement  of  the  geomagnetic  field  is  of  interest,  particularly  small  Gaus¬ 
sian  deviations  from  the  dipole  model.  Most  of  the  geographical  heterogeneity  has  been  smoothed 
out  at  this  height. 

2.  1000km  <  r  <  20Rg  :  This  includes  the  Van  Allen  radiation  belts.  The  field  is  generally  dipolar, 
being  approximately  300  nT  cl  5 Rb  and  10  nT  at  20/fs  (Re  is  the  radius  of  the  earth,  ca  0400 
km).*  Here,  deviations  due  to  intense  ionic  currents  produced  by  the  belts  arc  of  interest-  (Fig. 
2.7) 

3.  r>20RE:  Field  values  arc  less  well  known  in  this  region,  but  are  typically  5  nT.  Spatial  and 
temporal  variations  are  of  the  order  of  1  nT  and  are  due  to  phenomena  associated  with  the  solar 
wind. 

2.4  Magnetic  Units 

The  subject  of  electromagnetic  units  is  a  thorny  one  and  no  where  is  this  more  evident  than  for  magnetic 
units.  Historical  reasons  brought  about  the  plethora  of  units  which  has  led  to  considerable  confusion.  The 
reader  can  remove  much  of  the  mystery  concerning  magnetic  units  by  referring  to  the  excellent  discussion 
by  Jackson  (5).  Here  we  shall  only  list  the  units  which  are  normally  used  (Table  2.2).  Fortunately,  51 
units  are  now  considered  to  be  the  standard  for  all  but  certain  theoretical  physics  applications  and  we 
will  use  them  in  this  course.  SI  units  are  in  boldface  in  Table  2.2. 

*A«  a  crude  cttiwalc  of  the  maximum  angular  etianfe,  AS  ~  coa~>  (J||§^  )  —  l.T°. 

‘The  magnetic  tlorm  of  1SS9  produced  (A P/F)  ~  10%. 

*The  field  obey*  aa  e~*  la*. 
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2.5  Fields  of  Magnetic  Anomalies 

Deviations  of  the  measured  magnetic  field  due  to  sources  other  than  those  of  geomagnetic  origin  previously 
mentioned  are  usually  called  "anomaly  fields”  or  “magnetic  anomalies”.  Geological  anomalies,  due  to  ore 
bodies,  are  typically  in  the  range  of  100  -  1000  nT.  Virtually  all  magnetic  anomalies  are  due  to  magnetite 
concentrations  which  are  associated  with  the  ore  of  interest  (gold,  asbestos,  uranium-niobium). 

Magnetite,  being  ferrimagnetic,  has  a  magnetic  susceptibility  typically  between  1.2  and  13.  Ferro¬ 
magnetic  materials  have  susceptibilities  almost  always  greater  than  50  in  the  earth’s  field.  Man-made 
ferromagnetic  objects  have  much  smaller  dimensions  than  magnetite  formations  but  sensors  can  usually 
get  much  closer  to  them.  Thus,  on  the  whole,  the  fields  associated  with  ferromagnetic  objects  of  interest 
tend  to  be  of  the  same  order  of  magnitude  as  geological  anomalies. 

Typical  nongeological  objects  of  interest  may  include  armoured  vehicles  (peak  field  ~  10000  nT), 
submarines  (~  1  -  10  nT),  buried  mines  and  artillery  shells  (~  10  -  1000  nT),  archaeological  artifacts, 
parts  of  the  human  body  (<  10~*  nT)(Fig.  2.8). 

2.6  Methods  of  Detection 

Sensors  which  measure  magnetic  fields  are  called  magnetometers,  although  usually  the  term  is  reserved 
for  sensors  measuring  fields  less  than  the  typical  geomagnetic  field.  Generally,  there  are  two  types 
of  magnetometer  -  vector  sensors  and  total  field  sensors.  The  former,  which  include  fluxgate  sensors, 
optical  fiber  interferometers,  thin  film  sensors  and  Superconducting  Quantum  Interference  Devices  or 
SQUIDs  measure  a  single  component  of  the  field.  The  latter,  which  includes  all  magnetic  resonance 
magnetometers  (proton  precession,  optically  pumped  alkali  vapour),  measure  the  magnitude  of  the  field 
but  not  its  direction.  Fig.  2.8  shows  minimum  sensitivities  of  typical  magnetometers  and  clearly  several 
of  these  are  sufficiently  sensitive  for  our  purposes. 

Principles  of  operation  of  various  magnetometers  are  beyond  the  scope  of  these  lectures,  but  this 
information  may  be  found  in  the  bibliographic  references.  As  an  example,  the  operation  of  a  commonly 
used  vector  magnetometer,  the  fluxgate,  is  illustrated  in  Fig.  2.9. 

In  general,  if  the  anomaly  field  is  denoted  by  b  =  (bi,b},6j)and  the  earth's  field  is  6o  ss  (5<n,boi,bas). 
then  a  vector  magnetometer  would  measure 

^*=4^+6,  j=l,2or3  (1) 

whereas  a  total  field  magnetometer  would  measure 

h,  a:  j£  +  b|»  +  b*  +  ■  $)''*  (2) 

Generally  fco  >  b  and  a  Taylor  expansion  can  be  performed.  We  get 

i  [*>-  M’]  -  JgH  [v-M’J  +•••  <» 

7  =  t>0%  (4) 

Thus  to  first  order,  a  total  field  magnetometer  measures  the  magnitude  of  the  earth's  field  plus  the 
projection  of  the  magnitude  of  the  anomaly  field  along  the  direction  of  the  earth's  field  vector. 

The  chief  advantage  of  total  field  magnetometers  is  that  they  are  insensitive  to  small  changes  in  the 
orientation  of  the  tensor.  Vector  magnetometers  will  suffer  targe  baseline  fluctuations  when  rotated  by 
only  a  small  angle,  due  to  the  large  ambient  field  (Fig.  2.10). 

One  way  of  eliminating  such  baseline  fluctuations  is  to  use  two  sensors,  precisely  aligned  and  spaced 
far  enough  apart  that  different  anomaly  field  values  exist  at  the  different  sensors.  Such  an  instrument  is 
called  a  gradiometer  (although  a  true  gradient  is  measured  only  if  the  sensor  spacing  is  much  lets  than 
the  sensor-to-object  distance).  Alignment  of  sensors  is  very  critical  as  Fig.  2.10  shows. 

If  3  vector  magnetometers  are  situated  close  together  and  oriented  to  measure  3  orthogonal  components 
of  the  field,  they  can  be  used  as  a  total  field  magnetometer.  One  advantage  of  this  method  is  that  the 
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baseline  Held  6<y  may  be  subtracted  from  each  jth  component  by  extrapolating  from  an  anomaly  free 
region  and  thus  h|  may  be  measured  instead  of  |  bo  +  hjj.  Generally,  the  3  sensors  are  made  orthogonal 

only  to  ~  0.01°  and  software  compensation  is  applied  to  effectively  improve  the  alignment. 

Baseline  drift  due  to  temporal  fluctuations  of  the  earth’s  Held  may  be  eliminated  by  the  use  of  a 
“baseline”  or  reference  magnetometer.  Simultaneous  field  readings  are  obtained  both  for  the  reference 
magnetometer  and  another  which  is  used  to  measure  the  field  near  an  anomaly.  If  the  reference  magne¬ 
tometer  is  situated  far  from  any  anomalies,  its  field  reading  may  be  subtracted  from  the  reading  of  the 
second  magnetometer.  Spatial  uniformity  of  the  geomagnetic  field  (Section  2.2)  ensures  that  geomag¬ 
netic  fluctuations  will  be  eliminated  from  the  second  magnetometer  signal  by  such  a  method.  Baseline 
elimination  using  total  field  magnetometers  is  slightly  trickier  if  very  high  precision  is  required  due  to 
the  nonlinear  nature  of  Equation  2. 

High  accuracy  magnetic  measurements  must  be  performed  in  a  low  magnetic  noise  environment.  There 
are  4  ways  of  achieving  such  an  environment. 

1.  Low  metallic  content  laboratories!  These  are  situated  in  remote  locations  (at  least  40  km  from 
high-voltage  power  lines  and  a  few  km  from  other  buildings).  There  is  very  little  metal  present  and 
none  of  it  is  ferrous.  All  power  lines  are  shielded.  In  auch  locations,  field  gradients  as  low  as  0.1 
nT/m  are  obtainable.  The  schematic  for  a  low  metallic  content  laboratory  is  shown  in  Figures  2.11 
and  2.12. 

2.  Magnetic  shields!  These  usually  consist  of  layers  of  high  permeability  material  with  insulating 
gaps  between  them.  Gap  spacing  and  material  thickness  are  chosen  to  optimise  field  attenuation. 
A  properly  designed  shield  can  reduce  the  interior  field  by  as  much  as  IQ-®  over  the  external  field. 
Whole  rooms  as  well  as  small  cylindrical  volumes  can  be  shielded  in  this  manner. 

3.  Active  compensation;  A  field  free  region  may  be  obtained  by  measuring  three  orthogonal  com¬ 
ponents  of  the  external  field  using  a  magnetometer  and  applying  an  electrical  current  to  a  set  of  3 
Uelmholts  coils  with  mutually  orthogonal  sues  to  produce  a  field  opposite  to  the  external  field. 

4.  Gradiorueter  or  double  gradiometur  measurement!  If  it  is  desired  to  measure  the  very  weak 
field  of  a  source  in  close  proximity  to  the  sensor,  one  can  measure  the  derivative  or  second  derivative 
of  the  field  using  a  gradiometcr  or  double  gradiometer  respectively.  Uiomagnetic  measurements, 
such  as  measuring  the  field  of  the  heart  or  magnetite  in  lungs,  arc  typical  applications  of  this 
technique.  Roughly  speaking,  the  field  of  a  dipolar  source  falls  off  as  r~3,  r  being  the  souree- 
to-sensor  distance.  The  derivative  and  second  derivative  fall  off  as  r~*  and  t~%  respectively  and 
thus  unless  a  source  is  close  to  the  detector,  its  measured  gradient  or  double  gradient  will  be  quite 
small.  For  example,  if  a  source  at  1  cm  from  a  sensor  is  moved  to  10  m,  its  measured  field  is 
reduced  by  10"*,  its  gradient  by  1Q-U  and  its  double  gradient  by  10”**.  Biomagneiic  sources  have 
been  measured  using  double  gradiomelers  in  unshielded,  uncompensated  rooms  immediately  above 
subway  stations. 


2.7  Rotational  Properties  of  Magnetostatic  Measured  Quantities 

2.7.1  Rotationally  Invariant  Measurements 

The  most  commonly  measured  quantities  derived  from  the  magnetostatic  field  vector,  B,  arc  the  three 
components,  D<,  the  magnitude  of  the  field,  B,  the  gradient  tensor,  G  and  the  double  gradient  tensor, 
r.  The  i,j  component  of  the  gradient  tensor  is 


where  a  cartesian  system  with  position  given  by 

A?r^(A'»,  *,.*,) 


(6) 
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is  assumed  (superscript  T  denotes  the  transpose).  Likewise  the  double  gradient  is  defined  by 


r  -8G<1 

r«*  -  ixi 


m 


Because  measurement  platforms  are  seldom  stable,  it  is  very  important  to  know  how  these  quantities 
behave  under  a  rotation  of  the  coordinate  system. 

If  a  rotation  fron  an  unprimed  to  a  primed  coordinate  system  is  defined  by  an  (orthogonal)  rotation 
matrix,  A,  whose  elements  are  the  direction  cosines  for  the  transformation,  then  the  vector  B  transforms 
according  to  [9] 

B'  =  A  B  (8) 

G  is  a  second  rank  tensor  and  transforms  according  to 


G'  =  AGAt  =  AGA"1 


(») 


We  wish  to  determine  what  quantities  exist  that  are  invariant  under  a  rotation.  Clearly  the  magnitude 
B ,  being  a  scalar,  is  invariant.  Invariant  quantities  associated  with  the  gradient  tensor  become  apparent 
when  we  attempt  to  find  its  eigenvalues  and  eigenvectors.  To  do  this  one  must  set  det(G  -  AI)  =  0, 
where  A  is  an  eigenvalue  and  I  is  the  3  x  3  unit  matrix.  It  is  easy  to  show  that  det  (G  -  AI)  is  an  invariant 
quantity  for  any  A  and  thus 

dct(G'-AI)  =  det(G-AI)  (10) 

Direct  evaluation  leads  to 

det(G-AI)  =  D  +  AQ  +  A,T-A3  (11) 

where 

D  =  (GuGjjGjs  +  GisGmCji  +  GuGjiGjj) 

-  (GuGjsGjj  +  GiaGjjGss  +  GtsGjiGjj)  (12) 

<?=-^Tr  (GrG)  « 

(GijGjt  GuGji  +  GjjCjj)  -  (GijGh  -f  G\yGu  +  G«GM)  (13) 

T  =  Gu  4  Gjj  +  Gjj  (14) 

Inspection  reveals  that  D  is  the  determinant,  Q  is  the  negative  sum  of  principal  axis  cofactors  and  T  is 
the  trace  of  G. 

By  substituting  Equation  11  in  10  wc  get 


and  since  A  is  an  arbitrary  constant  we  have 


From  Maxwell's  equations  we  have 
But  Trace  (G)  =  V  •  B  and  hence 


which  implies  that 


{Q-‘Q,)A-X1(T-r)  =  0 

(15) 

D-D’ 

(16) 

Q  =  Q‘ 

1! 

*4 

< 

V-B- 0 

(17) 

T  =  0 

(18) 

Maxwell's  equations  also  state 

V  x  B  =  0 

(19) 

f  =  GJt  *  ^  3 

(20) 

e 


The  field  gradient  tensor  is  thus  symmetric  and  is  specified  by  only  5  independent  elements  (any  two  of 
the  diagonal  elements  and  either  the  upper  or  lower  off-diagonal  elements).  The  invariant  forms  D  and 
Q  can  be  reduced  to 

D  =  (Gu  -f  G33)  (G|j  —  G32GJ3)  +  2G12G13G23  —  (Gj2G33  +  GJ3G22)  (21) 

Q  =  Gj2  +  +  G33  +  G33  -f  G^  +  G23C33  (22) 

In  order  to  investigate  the  properties  of  D  and  Q  we  choose  a  rotation  that  diagonalises  G.  This  is 
called  a  “principal  axis  transformation”  and  the  coordinate  system  axes  are  called  the  “principal  axes”. 
We  denote  this  coordinate  system  by  a  double  prime  (''). 


/iOj-1)  0  o\ 

g"  =  g^3  0  -|(n  +  i)  0 

\  0  oi/ 


where  the  asymmetry  factor  q  is  given  by 


V  = 


(23) 


(24) 


and  G(\,  Gjjt  G'^  are  the  principal  axis  components  of  the  gradient  tensor.  Without  loss  of  generality 
we  choose 

!G«I  >  |G»I  >  |G'ul  (25) 

and  thus 

0<u<  1  (2ft) 

Substitution  of  Equations  23  and  24  into  21  and  22  yields 


i(l-,,»)G'*  (2?) 

I  (3  +  ^)6^  (28) 

Equations  26,27 ,28  yields 

0  <  £>  <  (29) 

SQ<  0%  (30) 

Although  there  would  appear  to  be  only  2  independent  quantities  nteded  to  construct  G"  in  this 
system,  there  are  indeed  S,  since  3  Euler  angles  arc  necessary  to  specify  the  rotation  relative  to  a  known 
space-fixed  system  (usually  the  detector  system).  Mote  specifically,  the  transformation  matrix  AM  be¬ 
tween  the  principal  axis  and  space-fixed  system  is  formed  by 


A"  =  (e‘i  e*j  <3)  (31) 

where  e*i,«Y e*j  are  the  normalised  eigenvectors  of  G. 

Both  D  and  Q  require  all  5  independent  components  ul  u  and  both  are  rotationally  invariant.  Q, 
having  only  quadratic  as  opposed  to  cubic  terms,  should  tend  itself  to  easier  analysis  when  attempting  to 
determine  the  source  characteristics  from  the  field.  Q  can  more  closely  approximate  the  largest  principal 
axis  component  of  G  because  Qx^  has  a  smaller  range  of  values  than  D1^,  for  a  given  value  of  G^3.T  All 
these  points  suggest  that  Q  should  be  preferred  to  D  when  attempting  to  analyse  gradient  data.  This 
indeed  turns  out  to  be  the  case  when  the  properties  of  D  and  Q  are  compared  for  the  specific  case  of  a 
static  magnetic  dipole  source. 

rNol»  Owl  0  <  sad  0.W6G",  <  Ql  <  G^. 
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2.7.2  RotationaUy  Invariant  Measurements  of  a  Dipole  Source 

Now  lets  examine  the  magnetic  field  of  a  static  magnetic  dipole.  In  practice,  the  magnetic  fields  associated 
with  ferrous  bodies  frequently  behave  as  though  they  were  due  to  a  static  magnetic  dipole  situated  at 
the  geometric  center  of  the  body.  Our  aim  is  to  try  to  relate  the  invariant  measurable  field  quantities  to 
the  source  characteristics  (position,  dipole  moment)  so  that  the  latter  can  be  deduced  from  a  series  of 
field  or  gradient  measurements.  We  employ  two  coordinate  systems;  a  space-fixed  ore  denote  by  capital 

letters  and  a  system  fixed  to  the  dipole,  denoted  by  small  letters,  which  we  call  the  “body-fixed”  system. 

„  t 

The  dipole  is  located  at  position  Ao  =  (Aoi,  Aoa,  Aos)  in  the  space-fixed  system  and  at  (0, 0, 0)  in  the 
body-fixed  system.  The  well  known  expression  for  the  dipole  field  b  at  position  r7  —  (z^zj, zj)  in  the 
body-fixed  system  due  to  a  static  dipole  moment  mT  =  (mtlmj,T7ij)  is  given  by 

*  =  ^r"3  (-*  +  ^  l'9™]  f)  (32) 

where  po  =  4r  x  10-T  H/m  is  the  permeability  of  free  space  and 


D  the  body-fixed  system,  the  jh  element  of  the  field  gradient  tensor  is  given  by 


= 


_  3bj 


dxi, 


Substitution  of  Equation  32  in  the  above  equation  yields 

3p<j 


m  =  tr 


Zjm»  +  z*mj  + 

where  is  the  Kronecker  delta,  or  in  tensor  notation, 


■(J 


~  “3“ 


3 


8  =  (fmr  +  iftr*  +  [frin]  [i  -  *->«*]) 


(33) 


(34) 


Note  that  Equation  33  has  built  into  it  the  fact  that  the  gradient  tensor  i*  symmetric  and  traceless. 

The  space-fixed  and  body-fixed  coordinate  systems  arc  connected  by  a  rotation  matrix,  A,  such  that 

=  A3?  (35) 

The  same  transformation  applies  to  ail  other  vectors  as  well.  One  can  also  define  1)  and  G  in  the 
space- fixed  system. 

/  «  .  .  \ 

(30) 
(37) 


On 


-£*"( 


11  j  if  4  +  /t*  Afj  + 


ftwsl) 

list  J/ 


where  Rr  =  (.Y  -  A0)T  and  RJ  as  Rr  !L  • 

Now  we  wish  to  evaluate  some  rotatioaally  invariant  measurable  quantities.  The  square  of  the  field 
can  in  principle  be  obtained  from  a  3-axis  total  field  magnetometer. 


ba  *  PS  =  (g) : 1  r“*  (mT*n  +  3  [mTr]  *  r‘») 


(38) 


To  evaluate  D  and  Q  it  is  easiest  to  choose  a  body-fixed  frame  in  which  »r»i  -  mj  =  0  and  m»  =  m. 
D  and  Q  are  then  evaluated  using  Equations  21,  22  and  33.  The  result  is  expressed  in  an  explicitly 
invariant  form. 


D  = .  (£)*,-  K.)  (|~of  +  [i'i]  1^1) 


(39) 


*  C«a  jou  bk  ttt  equation*  (at  (,  8,  g  ud  G  is  it«v  lui  S  =  A B  ud  g  s  AGArt 


a 


(40) 


Now  assume  a  coordinate  system  in  which  m2  =  0  and  in  which  the  *3  axis  is  parallel  with  the  X;< 
axis.  Further  assume  that  the  zj  axis  makes  an  angle  a  with  the  X2  axis  (Fig.  2.13).  Field  measurements 
are  made  in  the  X3  =  z  plane  orthogonal  to  the  Z3  axis.  This  is  a  geon  etry  common  to  many  magnetic 
source  detection  problems  such  as  the  detection  of  buried  artillery  shells  and  mines,  submarines,  etc.  The 
measurement  plane  in  most  cases  would  be  horizontal  and  the  *3  axis  would  be  vertical.  The  rotation 
matrix  becomes 

(cos  a  sin  a  0  \ 

-sin  a  cos  a  0  (41) 

0  0  1/ 

As  an  example  of  the  analysis  of  an  invariant  quantity  we  will  examine  Q.  A  number  of  properties 
of  Q  in  the  measurement  plane  can  be  examined,  such  as  zeros,  full  width  at  half  maximum  (FWHM), 
maxima,  minima,  etc.  We  first  study  the  position  and  value  of  the  extrema.  Extrema  of  Q  in  the 
measurement  plane  occur  for 


8Q_ 

8xj 


=  0  for  j  — 1,2  simultaneously 


In  the  present  coordinate  system  Q  (Equation  40)  can  be  rewritten 

^  4  +  03*1  +  <*4*1*3) 

9=00  (»!  +  »!  +  »!)’ 


(42) 


(43) 


where 


<*o 


=  (4?)  !  “i  =  3t7*i  +  ™3 


aj  =  ;  a$  =  m*  4  3mJ  ;  04  =  4r.»ims 

Substitution  of  Equation  43  into  42  and  extensive  manipulation  shows  that  extrema  with  respect 
occur  in  the  *3  =  z  plane  for 


*lm«a  — 


3r?ii  +  m\ 


/minis,  ,  ,\l/S  (m  , \l/s  3 

(—h  +  h)  +(-5-/1-/.)  -j’"1"5 


,»/*  3 


where 


and 


fi  =  27 4  Slrnfml  4  23mj 

/a  =  [c0mJJ  4  cim}°m|  4  4  c3mjmj  4  4  4  camJJ]l/J 

1  / 1836  729  \  ( 2754  6714\ 

e°  -  -  ;  cx  -  ^i?28  +  04  J  :  Cj  =  y  64  +  1728  J 

_  / 3843  13481 \  <  _  /2346  156Q6\ 

°3  ”  \  64  +  1728  j  ’’  “  V  64  +  1728  ) 

( 529  9996  \ 

Cs=  V  64  +  1728  J  ;  Cfl  = 


(44) 
to  *1,  aj 

(45) 

(46) 

(47) 

(48) 


2744 

1728 


(49) 


Substitution  of  *Sm««  into  the  equation  for  Q  (43)  yields  a  maximum  value  of  Q  of 


(mm  —  *  *<*0 


(flS  4"  ~  [<*4  +  flijy  ^1  4 


(50) 


where  at  is  the  quantity  in  [  ]  brackets  in  Equation  45.  The  subscript  umazn  has  been  used  above  to 
denote  extrema  points  since  it  can  be  straightforwardly  shown  the*  Q  >  0. 

Inspection  of  Equation  45  reveals  that  there  is  only  a  single  maximum  and  tl  at  its  position  is  on 
the  axis  defined  by  the  projection  of  the  horisontal  component  of  'he  d;pole  moment,  displaced  from  the 
origin  by  an  amount  that  is  a  fraction  of  z  (Fig.2.14).  The  fraction  is  a  furction  of  Tti/mj  and  varies 
from  0  at  mi  =  0  to  a  maximum  as  0.15z  at  mi  as  2.5ms  and  back  to  0  at  mj  =  u.  Thus  the  maximum 
is  directly  above  the  dipole  to  within  <  15%.  This  is  illustrated  in  Figs.  2.15-2.17  where  maps  of  the 
value  of  Q  in  a  plane  lm  above  a  dipole  (z  =  lm  plane)  are  shown.  The  Figures  correspond  to  different 
dipole  orientations.  Fig.  2.15  corresponds  to  a  dipole  with  mi  =  0,  —  1A  m*,  Fig.  2.16  to  a  dipole 

with  mi  =  1,  tnj  =  OA-m1  and  Fig.  2.17  to  a  dipole  with  mi  =  2,  ia$  —  lA-m1.  For  comparison, 
Figs.2. 18-2.20  respectively  show  the  maps  of  for  the  same  dipole  orientations  and  geometry.  Note 
that  the  map  for  Q  is  appreciably  more  narrow  than  that  for  63,  i.e,  the  field  strength  is  concentrated 
nearer  to  the  position  directly  over  the  dipole.  Furthermore,  note  that  there  is  a  single  maximum  for  Q 
which  is  almost  directly  over  the  dipole  whereas  there  are  maxima  and  minima  for  63  which  can  be  far 
removed  from  the  position  in  the  plane  immediately  above  the  dipole.  It  appears  that  using  Q  to  localise 
the  position  of  the  dipole  in  the  measurement  plane  would  require  much  less  signal  processing  than  63, 
provided  moderate  accuracy  is  required. 

To  determine  other  pertinent  quantities  such  as  the  moment  components  and  the  depth  is  not  so 
simple.  Since  the  extremum  solution  provides  3  independent  pieces  of  information  (why?),  at  least  three 
more  are  required  to  solve  for  the  3  position  components  and  3  moment  components  of  the  dipole. 
Additional  useful  information  could  include  the  FWHM  or  FWQM  (full  width  at  quarter  maximum). 
The  position  for  the  half  maximum  points  is  found  by  substituting  Qmam(2  (Equation  50)  in  Equation  43 
and  solving.  Unfortunately  solutions  must  be  obtained  using  a  root  extractor  algorithm  to  obtain  values 
of  *j,  aj.  The  solution  possesses  radial  symmetry  about  as  -  (0, 0, 0)  only  when  mi  =  0,8  For  other  dipole 
orientations  the  locus  of  the  half  maximum  is  not  circular  but  something  like  a  distorted  ellipse  which  is 
not  centered  about  aj  =  »j  =  0  in  the  plane  (Fig.2.21).  The  ratio  of  major  and  minor  axeal0(a/6)  of  the 
locus  is  a  function  of  mi  /ms  (Fig.2.22)  and  can  be  used  t<.  determine  the  latter  ratio.  The  length  of  the 
axes  are  a  function  of  x  and  can  be  used  to  estimate  the  depth  by  employing  a  calibration  curve  similar 
to  Fig.2.22.  Knowing  m\/m3  and  x  allows  us  to  determine  a  more  accurate  estimate  of  the  position  of 
*i  =  0  from  Fig.2.14.  Lastly,  the  magnitude  of  m  can  be  obtained  by  substituting  the  measured  value  of 
Qn*mi  the  estimated  values  of  mt/ma,  s,  xim«a,  and  using  =  0  in  Equation  43. 

A  similar  analysis  may  also  be  applied  to  the  determinant  D  of  the  gradient  tensor. 

An  analysis  has  been  applied  to  the  3-axii  component  of  b  and  has  led  to  a  successful  nonrccursive 
algorithm  for  estimating  the  position  and  dipole  moment  of  a  static  magnetic  dipole  to  within  a  few 
percent  [31].  The  method  is  much  lets  complicated  than  the  above  method  for  analysing  Q,  since  the 
only  measured  quantities  required  are  the  field  values  and  positions  of  the  maximum  and  minimum  in  the 
measurement  plane.  An  added  advantage  is  that  63  will  be  much  larger  than  the  field  differences  required 
to  measure  the  gradient  and  so  a  lest  sensitive  sensor  is  required.  It  should  be  pointed  out  that  although 
63  is  not  in  itself  a  rotationally  invariant  quantity,  in  low  to  mid-latitudes  it  is  reasonably  approximated 
by  (b,  ~  bo)  and  hence  can  be  measured  by  a  total  field  magnetometer  whose  readings  sure  insensitive  to 
sensor  orientation. 

A  comparison  of  some  of  the  properties  of  Q  and  63  in  the  measurement  plane  it  provided  in  Table 

2.3. 

There  ate  also  nonrccursive  tracking  algorithms  based  on  gradient  measurements.  These  require  an 
eight  element  tensor  array  to  measure  the  five  independent  components  of  the  gradient  plus  at  lcasi  one 
field  component  [29]. 

2.8  The  Magnetic  Field  of  Permeable  Objects 

The  problem  in  magrutostatic  ouasi-remote  sensing  is  to  determine  the  nature  of  a  magnetically  permeable 
object  from  its  anomaly  field,  in  order  to  be  able  to  do  this,  it  it  first  necessary  *.0  understand  how  to 

•  Kw  this  special  case,  the  local  of  point*  who**  value*  arc  h*lf  ths  maximum  it  circular  with  •  radius  tt  0.40004*. 

“■The  locus  it  not  an  ellipse  but  we  cut  iUll  define  the  “ma]o»  axil"  to  be  the  symmetry  axis  and  tbs  “minor  axis*  to 
be  the  ocis  of  largest  length  orthogonal  to  the  major  axis  la  the  plane. 
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calculate  the  field,  given  the  magnetic  characteristics  of  the  source.  A  complete  discussion  is  beyond  the 
scope  of  these  lectures,  althcvjh  e_  jd  accounts  are  found  in  [5],  [8]. 

In  magnetostatics,  the  basic  assumption  is  that  there  is  no  time  rate  of  change  of  charge  density  p 
anywhere  in  space.  Thus  the  equation  of  continuity  giving  the  current  density  J  becomes, 

V  •  J  =  0  (51) 

In  addition, 

VxH  =  f  (52) 

We  always  have 

V  •  B  =  0  (53) 

A  solution  of  Equations  52  and  53  that  satisifies  the  boundary  conditions  is  guaranteed  to  be  unique. 
Suppose  we  have  2  regions;  region  l  with  Si,  H\  and  permeability  pi  and  region  2  with  Bj,  //j,  pj.  Let 
n  be  the  unit  vector  directed  from  region  1  to  2  cn  the  surface  connecting  the  two  regions.  If  K  is  the 
idealised  surface  current  density,  then  the  boundary  condition  equations  are 

B\  •  n  =  Bj  •  h  (54) 

n  x  -  Hi)  =  K  (55) 

There  are  several  methods  of  solving  boundary-value  problems  in  magnetostatics.  They  are  broadly 
classed  as  follows: 

1.  Generally  applicable  method  •  vector  potential!  One  can  always  express  B  in  terms  of  a 
vector  potential. 

B  s  V  X  A  (56) 

We  need  a  constitutive  relationship 

B-B  (f?)  (57) 

but  unless  it  is  a  very  simple  one,  the  problem  is  intractable.  The  most  common  choice  is, 

B  =  ptf  (58) 

Substitution  in  Equation  52  and  choosing  the  Coulomb  gauge  •  A  =  oj  gives 

V*A  ~  pJ  (59) 

This  can  then  be  solved  subject  to  the  boundary  conditions. 

2.  Method  for  J  —  0  -  scalar  potential:  In  this  case 

H  =  -Vtfm  (00) 

where  <j>m  is  the  scalar  potential.  Again  we  need  t  constitutive  relationship  and  the  problem  is 
intractable  unless  Equation  57  is  simple.  If  Equation  58  is  valid  then 

=  0  (61) 

and  we  then  apply  boundary  conditions. 

3.  Bard  ferromagnots:  In  this  case  M  is  given  and  J  =  0.  It  can  be  shown  that 
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where  the  field  coot  'inates  are  unprimed  and  the  source  coordinates  are  primed.  It  is  assumed  that 
H  is  finite  within  a  volume  V1  bounded  by  surface  S'.  11 

Alternatively,  if  we  choose  to  use  the  vector  potential  A,  we  have12 


Note  that  for  aU  three  methods  the  problem  becomes  extremely  complicated  unless: 

1.  a  simple  B  ^7  j  relationship  exists  or 

2  a  simple  distribution  of  M  exists  and 

3.  the  geometry  is  very  simple. 

Only  a  few  special  cases  have  analytic  closed  form  solutions.  Among  these  aTe  the  uniformly  magnetised 
tphcre,  the  Homogeneous  ellipsoid  (which  includes  the  sphere  as  a  special  case)  in  a  uniform  external  field 
and  the  infinite  length  homogeneous  cylinder  in  a  uniform  external  field. 

Before  we  study  a  specific  boundary  value  problem,  a  few  words  should  be  said  about  demagnetisation, 
a  phenomenon  which  is  commmon  to  virtually  all  magnetostatic  problems.  Simply  put,  the  induced 
magnetisation  of  a  permeable  body  is  not  given  by  the  product  of  the  external  magnetic  field  and  the 
volume  susceptibility,  but  is  reduced  by  a  factor.  It  is  simple  to  show  how  this  arises. 

Su,  pose  a  body  with  volume  susceptibility  \  is  placed  in  a  magnetic  field  of  intensity  Hq.  Inside  the 
body  there  will  be  a  magnetic  field  Hi  while  outside  there  will  be  a  field  H  =  Ho  +  H,  where  H,  is  the 
secondary  field  due  to  the  magnetisation  M  of  the  body.  When  one  solves  the  boundary  value  problem, 
one  often  finds  a  relationship  of  the  form 

Hi  =  H0  -  dM  (64) 


where  d  is  u  scalar,  flow  jr, 

&  =  XJ?i  (65) 

so  that 

Htr  H0/(l  +  dX)  (66) 

Thus  the  appa*enl  susceptibility,  that  is,  the  quantity  which  •.  ?hen  multiplied  by  the  external  field  gives 
the  magnetisation  is 

v  -  „JL„ 

*  1  +  dx 

The  terra  dx  represents  the  reduction  of  the  fluid  itride  the  body  due  to  its  magnetisation. 

The  factor  d  is  eeUcd  the  “demagnetisation  factor"  and  ranges  from  0  for  needle-like  objects  magne¬ 
tised  along  their  long  axes  to  1  for  a  Ante  plate  magnetised  transversely.  As  an  example,  a  uniformly 
magnetised  sphere  has  a  =  1/3.  Thu*  demagnetisation  reduces  the  magnetic  moment  of  a  body  by  the 
ratio  1/  (1  +  dx)  (0  <  d  <  1)  and  is  a  function  of  the  orientation  of  tb-  body  with  respect  to  the  ambient 
magnetic  field.  Obviouriy,  if  x  is  small,  demagnetisation  is  negligible.  This  is  typically  so  for  diamagnetic 
and  paramagnetic  materials.  Demagnetisation  may  be  of  considerrblc  importance  for  ferromagnetic  and 
ferri magnetic,  materials. 

With  these  general  considerations  in  mind,  we  turn  our  attention  to  the  homogeneous  spheroid  in  a 
uniform  external  magnetic  field.  Although  this  has  been  solved  -xactly  [8],  we  shall  approach  the  problem 
difictciitly  by  performing  a  multipole  expansion  analysis.  This  will  allow  us  to  demonstrate  the  nature  of 
the  various  significant  muUipole  components  of  the  field.  In  c^-.tion,  the  method  is  applicable  to  objects 
whose  shapes  arc'  more  complex  and  which  do  not  have  closed  form  solutions.18 

1 1  Note  that  -  V*  •  Si  =  Pm  where  i>m  1«  the  effective  magnetic  volume  charge  density  and  A?  •  n*  =  <r«  where  A  is  the 
unit  vector  paralld  to  the  surface  dement  dr'  ard  <r„  U  the  eff«.  stive  magnetic  surface  charge  deadly. 

11  A?  K  n'  tj  the  effective  magnetic  turl'ace  current. 

!,Thc  field  inside  a  susceptible  volume  is  only  uniform  if  the  lulcrior  properties  are  homogeneous  and  isotropic  and  the 
boundary  surface  is  tecotul  order,  such  eo  in  the  case  for  an  ellipsoid.  It  Is  convenient  wh*n  analysing  commonly  encountered 
shapes  such  sls  sheets  and  prisms  to  assume  a  uniform  interred  field.  The  aasun  plion  will  fail  near  the  volume  boundary 
but  is  usually  a  reasonable  approximation  through  the  body  as  a  whole. 
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2.9  Example  of  Magnetic  Modelling  -  Multipole  Expansion  Method 

Assume  a  uniformly  permeable  spheroid  in  a  uniform  parallel  magnetostatic  field.14  Assume  a  cartesian 
coordinate  system  defined  by  orthogonal  unit  vectors  (z[,  z'Jt  with  origin  at  the  geometric  center 

of  the  spheroid.  The  symmetry  axis  is  the  z'a  axis  (Fig.2.23).  It  can  be  shown  (see  [8], p.257)  that  the 
induced  field  inside  the  spheroid  is  also  uniform  and  parallel.  Inside  the  spheiuid  we  have 

Ml=V^lBl-H !  (68) 

where  M'  is  the  magnetization,  B  is  the  magnetic  flux  density,  H  is  the  magnetic  field  and  the  subscript 
“1”  denotes  interior  quantities.  If  the  spheroid  is  uniform  and  isotropic  with  relative  permeability 
then 

Bi  =  HrfloHi 

M' =(/*,-!)£  (69) 

Thus  M[  is  uniform  and  parallel  too.  In  the  following,  we  drop  the  subscript  “1”  from  M[  since  we  will 
be  referring  to  magnetization  explicitly  only  inside  the  spheroid.  Now  the  scalar  potential  of  a  source 
with  magnetization  Af'  is  given  by  [6]  14 


<t>m  (r)  = 

5 •*(?)*’ 

(70) 

By  expansion  of  the  gradient  we  have, 

4 --V'-M'dv') 

.  t"  Jv,  r '•  ) 

(71) 

But  V'  •  M'  =  0  so 

II 

*l~ 

(72) 

Now  we  carry  out  a  standard  multipolc  expansion  about  the  origin.14 

*•«-  s  h  [£  (?)l 

Recall  that  the  scalar  potential  of  a  2*-pole  with  2*-polc  moment  ro^  is  given  by  (0) 

*<»■)  =  ^  —— -  ( -)  (74) 

9  4xn!  dz^dz'p...  \r"J  v  ' 

Thus  comparing  Equations  73  and  74  the  first  four  moments  of  the  source  distribution  are,  iu  component 
notation, 

mi°)=  /  Id'  di' 

Js< 

•*Il  U  anumed  that  (he  elliptoid  hae  bo  permanent  m»zneti*ation. 

"Prime  tupercript  denote*  *oure*  coordinate*,  unprimed  denote*  field  coordinate*,  double  prime  denote  mixed.  Source 
volume  U  V  and  board! aj  turface  It  S'.  It  i*  implied  that  r  i*  sol  intide  the  bod;.  Why? 

"Summation  convention  U  u*cd  in  Equation*  73  and  T4,  i.e,  when  indice*  are  repeated  la  the  tame  term,  (ummalion 
over  tbete  indice*  U  Implied.  For  example,  =  Yl/s  m-/Vi 
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m 


TO 


TO 


.(») 

'aPi 


l>=  [  z'aM'-df 
Js> 

= /*■■*« 

=  X 


M'  •  d«' 


M'  •  d#' 


/S' 

,(U 


(75) 


where  m(°)  is  the  scalar  monopole  moment,  mi^  is  the  vector  dipole  moment,  is  the  quadrupole 
moment  (a  rank  2  tensor)  and  is  the  octupole  moment  (a  rank  3  tensor). 

As  an  example  of  how  the  integrals  are  evaluated,  we  calculate  to(°). 


TO' 


(0) 


-L 


M'  ■  dt‘  = 


f  M[z[  •  d*'  +  [  MjZj  •  d#'  4-  /  MjZj  •  d#'  = 
ds1  ^S'  ds* 

Af(  [  z[  •  ds'  +  Afj  /  * i  •  dV  +  MJ  /  = 

^s'  7s'  Js* 

•tydv'  +  M^J  (v'-a^dt/ =0 


(76) 


The  second  lost  line  follows  from  the  uniformity  of  AT'  and  the  last  line  follows  from  the  divergence 
theorem.17 

In  a  manner  analogous  to  Equation  76,  we  can  show  that  the  a  component  of  the  magnetic  dipole 
moment  m(l)  is1*  19 

to^  =  A#;v"  a  =  1,2,3  (77) 

Similarly,  the  quadrupole  moment  tensor  can  be  shown  to  be  identically  sero.  This  doe#  require 
the  assumption  of  uniform  parallel  M'  as  well  as  the  assumptions  of  symmetry  about  the  x'a  axis  and 
mirror  symmetry  about  the  plane  defined  by  *j  =  0  (fore-aft  symmetry).  Thus 


«S=» 


p«i 


The  assumptions  leading  to  the  last  equation  also  imply  that  all  even  order  moments  are  identically  sero, 
i.e,  for  n  an  integer, 

in<2u>  =  0  (79) 

The  next  higher  order  contribution  to  is  due  to  the  octupole  moment  tensor  mi*!.  The  method  to 
calculate  the  moment  components  is  tedious  but  similar  to  the  method  used  to  calculate  the  previous 
three  moments.  Although  there  are  27  possible  elements  for  mj^7,  examination  of  Equation  75  reveals 
that  there  are  only  10  independent  elements.  The  assumption  of  a 3- axis  symmetry  further  reduces  the 
number  of  independent  nonzero  elements  to  6.  These  arc: 

m<’>  =  3m£>  =  3Af,'/u 

TOjjj  =  3to^jj  =  3Mj/n 
m&s  =  3Afj/jj 

mns  =  "•&  -  Af'/„ 


,T  We  have  proved  thi>  for  ihe  case  of  constant  Ai  but  it  i>  tme  for  any  A?’.  Con  you  prove  it? 

UTM»  equation  Is  true  even  if  A?1  is  not  conal&nl  throughout  the  body. 

is  independent  «f  the  choice  of  origin  iu  the  expansion  of  Equation  73  but  this  is  not  true  for  higher  order  momenta. 
Sec  [Si  p.130. 
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ro&i  =  M[hz 

=  1 (80) 

where 

hi  =  \  I  (r'J  "  *a )  dv>  (81) 

l  Jyt 

Isa=  f  *3 \dv'  (82) 

Jv 

m*]c0  =  mal«  =  mflac  (83) 

and 

=  °  if  (84) 

So  far  the  analysis  has  applied  to  any  axially  symmetric  body  with  fore-aft  symmetry.  Now  assume 
that  the  body  is  a  spheroid  and  that  the  diameter  of  the  spheroid  along  the  symmetry  axis  is  2oe  and 
the  diameter  of  the  largest  axis  orthogonal  to  the  symmetry  axis  is  2a.  It  is  simple  to  show  that 

V'  = 

o 

r  s 

A, =  jjea 

h,  =  (85) 

We  still  need  to  know  the  value  of  Af{,  Afj,  Afj.  Assume  that  the  external  magnetic  field  (usually  the 
earth’s)  in  the  absence  of  the  spheroid  is  =  (hotiioii&os)-  A  solution  of  the  boundary  value  problem 
((8],p.207),  assuming  the  permeability  of  the  spheroid  is  and  that  of  the  surrounding  medium  is 

Af;=^l^6oi  J  =  1,2,3  (80) 

Fj  -  (tbi  -  +  -  Mrs]/ (2#4»J)  (87) 

where  ^ 

Aj  =  ase  f  (»  +  aj)”1  (t  +  a*)  l(i  +  oV)  1/1  da  j=  1,2,3  (88) 

Jo 

with  ai  =  aj  =  a  and  aj  -  ac.  Integration  of  the  previous  equation  and  simplification  yields 

Ai  =  Aj  =  e  (e  +  £)  (eJ  -  l)  "l  (89) 

As  =s  -2e  (e“l  +  E)  (e*  —  l)”1  (90) 


£*l«  («-[«* -I]1)  («»-!)"* 


for  e  >  1  (prolate  spheroid), 


E  =  ^arctan  ^e  {l  -  e1}  -  x/2^  (l  -  cJ)  ^ 


for  e  <  1  (oblate  spheroid),  or 


Ai  =  Aj  =  Aj  =  -  (93) 

for  e  =  1  (sphere). 

Note  that  each  function  Fj  is  just  the  demagnetisation  factor  for  the  boundary  value  problem  correspond¬ 
ing  to  an  ambicut  magnetic  field  ityz).  The  functions  Fj  are  graphed  in  Figs.2.24-2.20. 
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Wc  now  wish  to  calculate  the  magnetostatic  field  at  point  r  from  the  origin,  due  to  the  field  induced 
in  the  spheroid.  We  assume  that  the  measurement  point  is  in  free  space.  If  we  keep  terms  up  to  and 
including  the  octupole  term  Equations  73  and  74  imply 


4>m  *  #(J)  +  #(8) 


or  in  terms  of  the  magnetic  field, 


“Po1  (fcW+hW)  (95) 

Following  extensive  manipulation  of  Equations  73,94,95,  we  have  in  component  notation  using  summation 
convention, 

bi3)  =  +  3j,_1  *«)  (®8) 

and 

*>i6)  =  ^r“5  (3ma^  ~  15r"2  +  ZBr-*xaxfixyxfm^fJ  (97) 

All  field  quantities  have  been  derived  in  a  body-fixed  coordinate  system  which  the  observer  does  not  a 
priori  know.  Wc  now  introduce  a  space-fixed  coordinate  system  with  arbitrary  origin,  whose  vectors  are 
indicated  by  upper  case  letters.  This  is  the  system  to  which  measurements  are  referenced.  It  is  assumed 

that  the  spheroid  center  is  located  at  R0  =  (Jfoii-Xos,  -X03). 

The  body-fixud  and  space-fixed  coordinate  systems  are  connected  by  a  set  of  Euler  angles  (#,  8 ,  #). 
The  angle  8  is  the  angle  between  the  spheroid  symmetry  axis  and  the  vertical  and  #  is  the  angle  between 
the  projection  of  the  symmetry  axis  on  the  horiiontal  plane  and  the  space-fixed  1-axis.  Because  of  the 
axial  symmetry  of  the  spheroid,  wc  cun  choose  #  =  0.  The  Euler  rotation  tensor  [9]  is  then  given  by10 

(cost?  cos  $  cos  Sain#  -sin#  \ 

-sin#  cos#  0  j  .  (98) 

sindcos#  sin 8 sin#  cost?  ) 

The  body-fixed  and  space-fixed  vectors  are  related  by 

60  «  AB0 

m(‘)  =  AAffr) 
b  =  AD 
M*>  =  ABM 

6(»)  1=  ABW  (90) 


Do  —  [Box  <  D01,  B(xi)  = 

(Bo  sin 80  cos #o,  Bosin0o*m#o,  B0cos#o)  (190) 

where  80  is  the  polar  angle  and  #0  is  the  atimuthal  angle  of  B 0  in  the  space-fixed  cartesian  system. 

Examination  of  Equations  77  through  98  will  reveal  that  given  the  tire  and  shape  of  the  spheroid  (a,  e), 
the  magnetic  material  properties  of  the  spheroid  and  the  surrounding  medium  (prtiPo)i  the  location  of 
the  geometric  center  of  the  spheroid  (,Voi ,  A'oi,  A‘0j)  and  the  orientation  of  the  spheroid's  symmetry  axis 

with  respect  to  the  space  fixed  system  (0,  #),  we  can  calculate  the  magnetic  field  j  due  to  the  presence 
of  the  spheroid  that  a  sensor  would  measure  at  a  point  in  space. 

An  example  of  the  success  of  such  a  model  is  illustrated  in  Figs.  2.27  and  2.28.  The  first  Figure  shows 
the  measured  signature  (solid  line)  and  theoretical  prediction  (dotted)  for  a  0.08m  radius  mild  steel  sphere 

S0Nole  that  AtX  e  I  where  I  it  the  3x3  identity  tentur. 
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whose  center  lies  directly  under  the  sensor  path  at  closest  approach.31  The  depth22  was  measured  to  be 
0.73m.  Deviation  of  the  model  from  theory  has  been  shown  to  be  due  to  remnant  magnetization.23  The 
second  Figure  shows  the  measured  signature  of  a  105mm  howitzer  shell  (solid  line)  and  the  theoretical 
signature  of  a  spheroid  of  similar  size  and  shape  (dotted  line)  to  illustrate  that  the  model  is  applicable 
to  realistic  situations  involving  axially  symmetric  objects  too.  Depth  was  measured  as  0.96m,  the  shell 
symmetry  axis  was  tilted  47°  from  vertical  and  its  projection  in  the  horizontal  plane  was  43°  from 
magnetic  north.  Spheroid  parameters  of  a  =  0.06m  and  e  =  3.5  gave  the  best  fit  to  the  data.  Measured 
values  for  the  shell  were  estimated  to  be  o  «  0.05m,  e  «  3.5.24  Remarkably,  even  though  the  model 
is  only  an  approximation,  the  good  agreement  was  maintained  as  the  shell  was  rotated  and  placed  at 
different  depths  without  having  to  change  the  model  parameters.  The  slight  difference  between  the  two 
curves  of  Fig.2.28  is  likely  due  deviations  of  the  shell  shape  from  that  of  the  model  as  well  as  remnant 
magnetisation. 

This  brings  up  a  very  important  point  about  magnetostatic  modelling.  Remnant  magnetization  can 
be  very  difficult  to  model  since  its  strength  and  direction  are  often  unknown  for  an  individual  object. 
Generally  it  is  ignored  and  this  turns  out  to  be  satisfactory  for  many  applications  in  which  the  object 
has  not  had  a  “colourful”  magnetic  history.  This  assumption  is  not  generally  true  for  geomagnetic 
applications.  In  some  cases,  such  as  the  measurement  of  magnetite  in  miner’s  lungs,  the  material  can  be 
subjected  to  a  known  large  magnetic  field  prior  to  measurement,  thus  assuring  that  remnant  magnetization 
dominates  induced  magnetization.  The  orientation  of  a  particle  with  respect  to  the  premagnetizing  field 
must  be  known  if  the  anomaly  field  of  the  particle  is  to  be  modelled  successfully.25 

Although  analytical  models  such  as  the  multipole  expansion  method  provide  insight  into  how  the 
fields  behave,  they  are  of  limited  use  in  many  practical  situations  involving  permeable  objects.  This  is 
particularly  so  when  one  has  to  account  for  hysteresis  curves  of  B  versus  H  in  a  material,  the  effects  of 
remnant  magnetization  or  object  boundaries  that  cannot  be  expressed  simply.  In  such  cases  numerical 
solutions  of  the  field  equations  are  necessary.  One  may  solve  either  differential  or  integral  equations, 
depending  on  the  problem  at  hand,  by  numerical  techniques  such  as  finite  difference  and  finite  element 
and  a  number  of  st  andard  methods  are  in  widespread  use.  Discussion  of  such  numerical  techniques  are 
beyond  the  scope  of  these  lectures,  but  [19]  is  a  good  reference  article. 

2,10  Inverse  Magnetostatic  Problems 

2.10.1  Introduction 

The  actual  problem  with  which  we  arc  faced  in  magnetostatic  remote  sensing  is  the  inverse  of  that  in 
the  last  section.  The  problem  may  be  posed  as  follows:  “Given  a  set  of  magnetic  field  or  field  gradient 
measurements  referenced  to  a  space-fixed  system,  determine  the  position,  sise,  shape  and  material  prop¬ 
erties  of  the  object."  Problems  in  which  field  measurements  are  used  to  infer  properties  of  the  source  of 
the  field  are  called  “inverse  problems". 

To  study  the  properties  of  an  inverse  problem,  we  shall  use  the  the  uniform  spheroid  as  an  example. 
We  assume  that  the  accuracy  of  measurements  is  such  that  moments  of  higher  order  than  octupole  may 
be  neglected.  Note  that  the  problem  is  greatly  simplified  compared  to  the  case  of  a  body  of  arbitrary 
shape.  For  the  spheroid  there  are  9  unknown  parameters;  2  sise/shape  parameters  (a,e),  2  magnetic 
material  properties  (p*tifM)i  2  angles  which  determine  the  orientation  of  the  symmetry  axis  (0,  <£),  and 
the  3  position  eoordinete*  of  the  spheroid  center  (Xoi,  Agj>  -Xos)-  For  a  body  of  arbitrary  shape,  there  are 
22  independent  parameters;  3  components  of  m(l),  6  independent  components  of  m^2),  10  independent 
components  of  and  the  3  position  components  of  the  center  of  the  object.26  Nevertheless,  even  the 

21  Senior  travelled  S  to  N  la  straight  tine  in  the  horizontal  plane.  Experimental  tel  up  it  thown  in  Fig.S.ll. 

22  Distance  between  center  of  the  tcator'i  active  volume  sad  the  geometric  center  of  the  object  at  the  point  of  doeetl 
approach  of  the  tensor. 

21  How  would  you  go  about  showing  thilT 

2*  We  also  have  p.,  =  100,  ss  1,  #o  «=  1T°,  a  180°.  The  ateel  portion  of  the  artillery  ihell  it  by  no  meant  a 
tpheroid.  It  it  axially  symmetric  but  it  it  hollow  and  hat  an  angular  shape,  somewhat  like  a  cross  between  a  truncated 
apheroid  and  a  cylinder.  The  eatlmates  of  a  and  a  for  the  ihell  are  thus  only  approximate  quantities  based  on  measuring 
the  overall  length  and  width  of  the  shell. 

2*Typically  the  particles  are  modelled  as  spheres  or  randomly  oriented  spheroids. 

**  Where  is  the  Information  concerning  the  site, shape,  magnetic  properties  and  orientation  of  the  body  contained? 


17 


uniform  spheroid  inverse  problem  turns  out  to  be  quite  difficult  to  solve. 

In  practical  situations,  the  spheroid  problem  can  be  further  simplified.  For  ambient  field  values 
typical  of  those  on  the  earth’s  surface,  we  usually  have  300  <  fit  1  <  1000.  Typical  rocks  have  values 
of  Hrj  between  1.0  and  1.1.J7  For  this  range  of  and  fit j,  Figs.2.24-2.26  reveal  thax  the  Fj  functions 
which  determine  the  components  of  the  magnetisation  of  the  spheroid  are  insensitive  to  fiti  and  /t,j. 
Thus  fit i  may  be  fixed  at  say,  500  and  fit 2  =  1  without  much  error.  There  are  now  only  7  unknown 
parametersja,  e,  8,  <j>,  Jfoi,  Jfoii  Aos- 

It  is  immediately  obvious  that  information  is  necessary  from  multipoles  of  order  higher  than  dipole  if 
the  problem  is  to  be  solved  uniquely.  The  space-  fixed  components  of  B(J)  are  a  function  of  6  parameters; 
the  three  components  of  Ao  and  the  three  components  of  Aff1).  Thus  there  is  a  built  in  degeneracy  that 
keeps  us  from  uniquely  determining  the  7  spheroid  parameters  from  B(J).58  Note  that  for  the  special 
case  of  a  sphere,  there  are  only  four  parameters  and  dipole  field  information  is  sufficient  to  solve  the 
problem.  Indeed  it  can  be  shown  that  the  homogeneous  sphere  in  a  uniform  external  magnetic  field  has 
a  pure  dipole  anomaly  field.  This  means  that  if  only  dipole  field  measurements  are  available,  we  cannot 
distinguish  between  a  compact  orientabie  body  of  arbitrary  shape  and  a  sphere.  We  will  later  see  that 
in  spite  of  this  apparent  deficiency,  dipole  measurements  can  still  be  used  to  identify  spheroids. 

There  are  two  general  techniques  that  may  be  used  to  estimate  the  spheroid  parameters  from  mea¬ 
surements  of  B.  The  two  are  “model  fitting”  and  “pattern  recognition”.  Although  our  problem  is  quite 
specialised,  these  methods  are  generally  applicable  to  all  magnetostatic  inverse  problems.18 

2.10.2  Solution  by  Model  Fitting 

Model  fitting  involves  devising  a  mathematical  model  to  describe  the  secondary  magnetic  field  as  a 
function  of  source  parameters  and  then  performing  maximum  likelihood  estimation  (MLE)  to  determine 
the  parameter  values  that  best  fit  the  measurements.  Most  geometries  do  not  have  simple  analytical 
models  and  are  thus  not  ideally  suited  to  this  method.  It  is  possible  to  use  a  numerical  model  in  the 
MLE  procedure  in  place  of  an  analytic  field  equation,  but  this  mokes  the  method  very  computationally 
intensive. 

For  the  spheroid,  we  do  have  an  analytical  equation.  Unfortunately,  in  this  case  as  in  most,  the 
field  equation  is  nonlinear  in  the  source  parameters  which  means  that  a  nonlinear  estimation  technique 
must  be  used.  Such  techniques  arc  recursive  and  usually  require  substantial  computer  time  if  there  ate 
many  data  points  and/or  parameters.  Also,  the  initial  values  of  the  parameters  must  be  guessed  and 
a  technique  which  is  insensitive  to  the  initial  guesses  is  necessary.  If  the  uncertainties  in  the  data  are 
Gaussianly  distributed,  least  squares  fitting  is  an  MLE  method  and  a  particularly  good  algorithm  to 
use  is  Marquardt's  method  [30].  The  technique  is  reasonably  fast  and  produces  accurate  estimates  of 
parameters  if  the  model  is  correct40  Its  main  strengths  are  that  its  performance  (i.e,  speed,  accuracy 
of  estimates,  etc.)  is  very  insensitive  to  the  initial  parameter  guesses  aud  it  is  reasonably  robust.  All 
fitting  methods  tend  to  be  mote  prone  to  oscillations,  instabilities,  and  nonconvergence  as  the  number 
of  parameters  increases.  This  may  be  partly  corrected  by  increasing  the  number  of  data  points  but  this 
only  helps  if  the  points  are  “well  distributed"  and  then  only  up  to  a  certain  number  (typically  ~  a  few 
hundred).  Also,  the  more  data  points  that  arc  employed,  the  longer  the  time  required  for  convergence 
of  the  algorithm.  In  spite  of  all  these  pitfalls,  spheroid  parameters  have  been  successfully  estimated 
using  the  “dipole+octupole"  model,  For  typical  laboratory  measurements  (position  errors  ~  0.01m,  field 
errors  ~  l.OnT),  parameters  have  been  estimated  to  an  accuracy  of  a  few  percent  [30],  Unfortunately, 
convergence  problems  and  the  formation  of  fnlsc  solution  sets  due  to  local  minima  in  xS  space  also  occur 
occasionally  and  must  be  circumvented  by  heuristic  checks  in  the  computer  code. 

"Magnetite  it  an  exception  with  2  <  pr  j  <  1  *  approximately. 

91  Hi3)  cui  be  measured  by  measuring  !I  at  distances  far  enough  from  the  object  to  that  id*)  Cj  0  but  near  enough  to 
that  i!  can  sttU  be  measured  with  sufficient  accuracy. 

"Specific  example*  may  be  found  in  (<]. 

3oThere  are  other  algorithm!  that  arc  fatter  and  more  accurate  but  not  both  together.  These  attribute!  lend  to  oppoae 
one  another  to  that  Marquardt'i  algorithm  it  a  good  compromise. 
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2.10.3  Solution  by  Pattern  Recognition 

Pattern  recognition  involves  comparing  characteristics  of  a  set  of  magnetic  data  from  an  unknown  object 
with  that  horn  a  known  object  to  determine  if  the  two  objects  are  the  same.  A  full  discussion  of 
pattern  recognition  is  beyond  the  scope  of  these  lectures,  but  a  number  of  good  introductory  texts  are 
available.  An  excerpt  from  one  of  them  (Tou  and  Gonzalez  [13])  is  included  with  the  notes  to  provide 
some  introductory  terminology.  Geophysics  routinely  makes  use  of  “characteristic  curves”  to  determine 
whether  the  source  of  a  magnetic  anomaly  is  a  dyke,  plate,  dipole  doublet  or  some  other  simple  body. 
The  chief  problem  is  that  generally  magnetic  profiles81  for  a  single  object  change  with  the  choice  of 
sensor  path  and  object  depth  and  orientation.  Thus  huge  amounts  of  data  must  must  be  stored  even 
if  the  number  of  possible  objects  is  small.  Clearly  some  form  of  data  reduction  or  compression  (called 
“feature  extraction")  is  essential  to  make  pattern  recognition  feasible  for  magnetic  identification.  One 
plan  is  to  oversample  each  profile,  then  perform  a  Fast  Fourier  Transform  and  discard  higher  frequency 
information  to  form  the  feature  vector  for  the  profile  8i.  For  certain  geometries,  heuristic  extrapolation  or 
interpolation  techniques88  can  be  employed  to  minimise  the  number  of  different  depths  at  which  profiles 
must  be  obtained.  In  spite  of  these  and  other  methods  of  data  compression,  profile  matching  methods 
usually  involve  huge  libraries  of  characteristic  curves  or  feature  vectors. 

2.10.4  Pattern  Recognition  for  Compact  Axially  Symmetric  Objects 

In  a  number  of  realistic  situations  the  sensor-to-object  distance  is  sufficiently  large  that  multipole  fields 
of  order  higher  than  a  dipole  may  be  neglected.  As  previously  mentioned,  this  means  that  an  orientablc 
compact  object  of  arbitrary  shape  cannot  be  uniquely  identified  by  measurements  of  the  field.  In  practice, 
however,  the  number  of  possible  object  shapes  and  sizes  for  a  particular  problem  is  usually  finite  and 
small.  Uniqueness  will  then  be  ensured  if  the  fields  associated  with  only  those  objects  are  uniquely 
different.  Under  these  conditions,  location  and  identity  of  a  compact  object  can  be  reliably  determined 
using  an  algorithm  that  is  fast  and  requires  only  a  amail  amount  of  storage  for  the  feature  vectors  of 
an  object.  We  will  present  the  method  for  an  axially  symmetric  object  using  a  spheroid  ns  an  example. 
The  method  can  be  generalized  to  an  arbitrary  compact  body  but  it  is  a  tedious  and  not  particularly 
illuminating  exercise. 

The  first  step  is  to  determine  the  location  and  dipole  moment  components  of  the  object  relative  to  a 
space-fixed  coordinate  system.  This  may  be  done  by  any  one  of  a  number  of  published  algorithms  using 
field  [30],  [31]  or  gradient  [29]  information  or  by  a  method  such  as  that  outlined  for  Q  in  Section  2.7.2. 
In  performing  the  measurements,  care  must  be  taken  to  ensure  that  the  sensor  is  far  enough  from  the 
object  that  higher  order  multipole  fields  are  negligable. 

The  space-fixed  components  of  the  dipole  moment  induced  in  a  homogeneous  permeable  axially  sym¬ 
metric  compact  object  by  a  uniform  magnetostatic  external  field  will  vary  with  the  orientation  of  the 
object.  Since  two  angles  define  the  orientation,  the  locus  of  all  possible  dipole  moments  for  such  an 
object  is  actually  a  two  dimensional  surface  in  a  three  dimensional  space.  This  can  readily  be  seen  for 
the  spheroid  problem  of  Section  2.9.  Using  Equations  77  and  86,  the  magnetic  moment  induced  in  the 
spheroid  due  the  ambient  magnetic  field  can  be  rewritten 

m  =  ^‘VF£  (101) 

where 

/  Fi  0  0  \ 

F  =  0  Fx  0  (102) 

\  0  0  F»  j 

and  V  s  V'.  By  substituting  Equations  98  and  99  in  101  and  following  extensive  simplification,  we 
find  that  the  space-fixed  magnetic  dipole  moment  vector  can  be  expressed  explicitly  in  terms  of  only 
space-fixed  quantities: 

11 A  profile  W  a  Ml  of  field  mcuunuxoll  taken  along  a  one  dimensional  curve,  usually  a  straight  line,  In  a  plane  above 
the  object. 

** This  relies  on  the  fact  that  for  most  magnetostatic  problems,  the  profiles  are  slowly  varying  spatial  functions. 

u  usually  based  on  scaling  the  spatial  extent  of  the  profile  as  a  function  of  depth 
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M  =  Mi1)  = 


Mi 

M2 

M3 


(103) 


({.Pi  +  [F3  -  Fi]  sin1  8  cos2  +  {[P3  —  Fi]  cos  8  sin  8  cos  $}.Bos  \ 

{[F3  -  Fusin'*  8cos<f>sin<l)}Bo\  +  {[Fz~  Fi]cos0 sm8sin4>}Boa  j  • 

{[Fa  -  Fi]  cos  8  sin  6  cos  <p}Bo\  +  {Fi  +  [Fs  -  Fi]  cos1 0}Bos  / 

Examination  of  Equation  103  reveals  that  for  a  given  ambient  magnetic  field  and  surrounding  medium, 
M  is  a  function  of  the  spheroid  shape,  size  and  material  (a,  e,  /i,i)  and  two  continuous  orientation  pa¬ 
rameters  Also,  because  of  symmetry,  unique  values  of  M  occur  only  for  8  =  <j>  =  0,  0  <  9  <  r/2 

when  0  <  <j>  <  2*  and  8  =  t/ 2  when  0  <  <f>  <  it.  Two  magnetic  dipole  moment  surfaces  corresponding  to 
two  spheroids  of  different  shape  and  size  are  shown  in  Fig. 2. 29. 

Our  plan  is  to  use  the  dipole  moment  vector  as  the  feature  vector  to  classify  the  object.  Because  the 
feature  vectors  form  a  surface  or  manifold  in  the  feature  space  (M  space)  we  must  find  a  classifier  that 
can  determine  to  what  surface  a  test  vector  belongs.  The  plan  is  to  find  a  2-D  analogue  to  the  Nearest 
Mean  Vector  classifier,  which  determines  the  distance  from  a  test  vector  to  a  point. 

A  number  of  restrictions  and  assumptions  will  be  imposed,  but  these  are  all  applicable  to  the  problem 
at  hand.  We  will  consider  a  prototype  which  is  a  manifold  of  M  dimensions  and  a  function  of  M 
independent  parameters.  The  feature  space  is  of  dimension  N  >  M.  It  is  assumed  that  the  manifold  is 
well-behaved,  ie.,  has  no  singularities  or  discontinuities  on  the  portions  under  consideration.  The  manifold 
is  assumed  to  be  finite  in  extent,  but  may  be  open  or  closed.  It  is  assumed  that  a  given  set  of  values 
of  the  M  independent  parameters  maps  into  no  more  than  one  point  on  the  manifold.  However,  one 
point  may  correspond  to  several  values  of  a  particular  independent  parameter.  (The  manifold  is  said  to 
be  “degenerate"  in  that  parameter  at  that  point).  It  ie  convenient  to  use  Dirac  bra-ket  notation  applied 
to  real  vectors.  Thus,  the  row  vector  (asi, xj...,x„)  is  written  as  (*[.  The  column  vector  with  the  same 
elements  is  denoted  by  |z)  and  the  scalar  product  of  (a|  and  |y)  is  written  (x|  y). 

We  will  specifically  be  considering  a  pattern  class  whose  prototype  consists  of  an  M  =  2  dimen¬ 
sional  manifold  in  an  N  =  3  dimensional  feature  space.  The  prototype  is  a  function  of  two  continuous 
parameters,  8  and  <p.  This  is  illustrated  in  Figure  2.30. 

The  prototype  for  class  i  and  given  8,<f>  is  a  point  defined  by  the  head  of  a  vector  denoted  |mj  (0, $)). 
The  prototype  is  approximated  by  a  finite  number  of  “unit  cells".  For  an  M  —  2  dimensional  manifold, 
the  unit  cells  are  truncated  hyperplanes,  or  hypertriangles  (triangles  for  N  -  3  as  in  Figure  2.30), 
which  connect  sampled  points  |m,-(^*)  on  the  hyperplanc.  The  subscripts  j,  k  indicate  that  the  prototype 
feature  vector  is  evaluated  at  discrete  values  of  the  parameter*  8  —  8j  and  4>  —  4>ik-  Where  on  the 
manifold  to  choose  the  in  order  to  optimally  approximate  the  manifold  is  dependent  on  its 

structure  and  will  not  be  addressed  further,  Consider  the  region  of  the  manifold  for  which  0j  <  6  <  0j+i 
and  tfk  <  4>  <  +  This  region  of  the  manifold  may  be  approximated  by  two  hypcrtrionglcs.  One 

hypertriangle  passes  through  ,  |w«sj,s+i)  and  is  bounded  by 

=  Kj.»+»)  -  irosj.*) 


and  Kyi4)  -  K^x). 

The  other  hyperlrianglc  passes  through  jri^j+j,*) ,  +  +  and  is  bounded  by 

K»)  :=  Kj  +  m)  “  Iwxj  +  i.i  +  i) 

\v(j.k)  T  Kj,*n>  ~  Hj+m+i) 


and 


(104) 


(105) 
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By  allowing  k  to  span  the  values  corresponding  to  the  upper  and  lower  limits  of  6j  and  <ph  respectively, 
the  entire  prototype  manifold  for  class  »  may  be  approximated.  In  what  follows  the  subscripts  »,  j,  k  will 
be  dropped  where  it  is  unambiguous  to  do  so. 

To  find  the  distance  to  the  hypertriangle  bounded  by  the  vectors  |ti) ,  |u), 

|«)  -  \v)  we  first  form  an  orthonormal  basis.  The  orthonormal  basis  |a) ,  \fi)  is  constructed  using  the 
Gram-Schmidt  method  [15].  It  is  assumed  that  |tt)  and  \v)  are  not  colinear.  The  basis  vectors  are 

|o)  =  |tt)/(u|u)i  (106) 


Let  |  a)  be  a  test  vector  and  define 


IWJ,»)  =  I*)  ~  \”kj,h)  ■ 


(108) 


If  we  further  define 


ViJ.h)  *°  be  *ke  projection  of  jjftj,»)  onto  the  subspace  spanned  by  the  basis  and 
\dij,h)  to  be  the  vector  which  is  normal  to  the  subspace  and  which  passes  through  }»),  then 

M)  =  Iv)  “  l^'5)  (109) 


and 


V<*>)  =  (vl  1°)  +  <vl  0)  1/3)  — 


M  _*>.  M  ,  ((vl v >  ~  (vl «)  W  •)  /  M  *))  /. .  («|  v)  \ 

(»l «)  +  ^(»|  »)-(«|  «)’/(«| »))  lw'5Riu,J' 

(110) 

After  substantial  simplification,  we  find 

=p|a)+tf|w) 

(111) 

where 

_  _  (Vl  *)  («l  *}  -  (y|  u)  (u|  v ) 

W  v)  («|  1i)  -  (u|  V )* 

(112) 

and 

__  (vl «)  -?H  v) 
p  («l«) 

(113) 

Note  that  the  denominator  in  Equation  112  it  sero  iff.  |u)  and  |v)  are  colinear. 

If  |a)  is  a  sample  from  the  class  i  corresponding  to  the  region  of  the  prototype  manifold  bounded 

bY  lu)  i  lv)  1 1**)  —  |v),  then  estimates,  9  and  of  the  continuous  parameters  associated  with  |a)  may  be 
obtained  from 


$  =  9j+  (0i+v  -  9j)  (114) 

^  =  #*  +  Pu,»(^*+i  -  &).  (US) 

The  minimum  distance,  <Lj,s,  from  the  test  vector  to  the  hypertriangle  is  approximated  by 

4<ja)  *  •  (116) 
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The  previous  equations  calculate  the  minimum  distance  to  the  infinite  hyperplane  passing  through  the 
hypertriangle.  The  region  of  the  manifold  under  consideration  has  been  assumed  to  be  approximated  by 
the  hyperplane  only  within  the  boundaries  of  the  unit  cell  (hypertriangle).  Thus,  the  previous  distance 
equations  are  only  valid  if  |j^)  lies  within  the  boundaries  of  the  hypertriangle.  It  can  be  easily  shown 
that  this  is  true  provided  p  >  0,  q  >  0  and  0  <  p  +  g  <  1.  If  any  of  these  conditions  are  not  satisfied, 
then  is  replaced  by  the  minimum  distance  from  the  test  point  to  the  vectors  which  bound  the 
hypertriangle,  using  the  method  described  above  for  approximating  the  minimum  distance  to  a  one 
dimensional  manifold.  For  speed  of  implementation,  q  is  calculated  first  and  based  on  its  value,  the 
appropriate  calculation  of  is  carried  out. 

By  changing  all  quantities  to  primed  quantities,  Equations  106  to  116  may  also  be  used  to  find  the  min¬ 
imum  distance  to  the  hypertriangle  bounded  by  the  vectors  |u') ,  jv') ,  |t')  -  |u')  provided  Equations  108, 
114  and  115  are  changed  to 


I yij.k)  =  1*)  -  Hj+i.fc+l)  (117) 

0  =  *iH  +  «Us(0i-0i+i)  (118) 

4>  =  <t>k+\  +  p(j>  (0s  -  &+i).  (110) 

The  minimum  distance,  d\,  tom  the  test  vector  to  the  manifold  is  then  approximated  by 

di  =  min{dij,t,d{ii>}.  (120) 

The  test  vector  is  assigned  to  the  class  i  for  which  is  a  minimum. 

The  prototype  manifold  for  certain  values  of  one  of  the  continuous  parameters,  say  6t ,  may  be  indepen¬ 
dent  of  the  other  continuous  parameter.  In  this  case,  the  manifold  between  <$t>  d(+i  >  and  may  be 
approximated  by  a  single  primed  hypertriangle  instead  of  both  a  primed  and  an  unprimed  hypertriangle. 

As  an  example  of  the  abilities  of  the  classification  method,  noise-free  magnetic  moments  of  six  different 
spheroids  (Table  2.4)  were  generated  at  15°  increments  of  the  orientation  angles  8,  $  to  produce  a  design 
set.  Test  sets,  consisting  of  magnetic  moment  feature  vectors  generated  at  5°  increments  were  generated 
with  different  %  noise  levels  added  to  them.  Table  2.5  shows  the  probability  of  misclassification  for  the 
classifier  just  described  (Continuous  Parameter  Classifier).  The  continuous  parameter  classifier  theory 
and  comparison  with  other  classifiers  is  discussed  more  fully  in  (32) . 

3  Electrostatic  Methods 

3.1  Introduction 

Of  all  the  electromagnetic  properties,  the  conductivity  («r)  has  the  widest  range  of  variation.  Whereas 
under  normal  conditions  the  magnetic  permeability  (p)  varies  from  ~  p$  to  lOOOpo  and  the  permittivity 
(«)  varies  from  ~  «o  to  80<q,  <r  can  span  20  or  so  decades  (sec  Fig.3.1,  Table  3.1). 

Two  effects  complicate  the  electrostatic  process  in  addition  to  simple  ohmic  conduction.  First,  po¬ 
tentials  can  develop  in  the  medium.  These  are  caused  by 

1.  a  difference  in  chemical  potentials  of  minerals  at  the  interface  between  two  minerals  or 

2.  gradients  in  solute  concentrations  in  interstitial  water  or 

3.  Quid  motion  in  porous  materials. 

These  are  steady  state  effects  which  are  generally  referred  to  as  “spontaneous  polarisation*.  Second, 
charge  may  accumulate  at  the  interface  between  certain  minerals  due  to  the  flow  of  current  from  an 
external  source.  This  is  c&iied  “induced  polarisation*.  The  difference  between  the  two  is  that  the  former 
may  involve  current  flow  in  the  absence  of  external  voltage. 
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The  “d.c.  resistivity”  or  “electrical  impedance”  method  is  the  simplest  electromagnetic  measurement 
method  there  is,  at  least  in  principle.  Electrodes  are  insetted  in  the  bx.dy  to  be  measured.  A  potential 
across  them  establishes  a  current  field  in  the  body  and  the  finite  conductivity  of  the  body  produces  an 
electric  potential  which  is  measured  at  the  body  surface  by  two  electrodes. 

In  geophysics,  the  body  is  the  earth  and  the  field  is  perturbed  by  the  presence  of  subsurface  zones  of 
conductivity  which  differ  from  the  otherwise  presumed  homogeneous  conductivity  of  the  volume.  There 
are  two  types  of  exploration.  Vertical  exploration  involves  delecting  layered  structures  and  is  generally 
done  by  making  measurements  as  the  electrode  array  gradually  increases  in  horizontal  spacing  about  a 
fixtd  point.  Horisontal  exploration  is  used  to  look  for  horizontal  anomalies  such  as  ore  bodies  and  involves 
moving  an  array  of  electrodes  of  fixed  spacing  horizon»«Jly  along  the  ground. 

Electrical  impedance  methods  are  used  in  medicine  to  measure  certain  global  or  bulk  cardiac  parame¬ 
ters  and  intrathoracic  fluid  volumes  and  some  attempts  have  been  made  at  imaging  soft  tissue  and  bones. 
Such  measurements  are  possible  because  of  the  large  conductivity  contrasts  in  the  body.®4 

There  are  many  types  of  exploration  array  systems  and  we  shall  only  mention  a  few.  For  vertical 
exploration,  the  most  popular  is  the  Wenner  array  (Fig.3.2).  The  voltage  for  a  given  input  current  is 
measured  and  the  array  is  expanded  about  the  center  so  that  deeper  layers  have  more  effect  on  the 
potentials.  We  shall  return  to  this.  For  horizontal  exploration,  the  current  electrodes  are  often  fixed  at 
large  distances  and  the  potential  contours  are  mapped,  i.e.,(Fig.3.3) 

There  are  some  problems  associated  with  taking  resistivity  measurements.  One  is  contact  potential, 
the  potential  between  electrodes  in  the  absence  of  source  current.  It  is  due  to  electrochemical  emfs  on 
electrode  surfaces  or  material  interfaces  in  the  host  body.  One  solution  is  to  measure  potentials  with 
the  current  in  first  one  direction  and  then  another.  Thir  may  be  done  using  low  frequency  AC  current 
and  measuring  only  the  AC  component  of  the  voltage,  but  the  frequency  must  be  low  enough  so  that 
potentials  induced  by  the  magnetic  field  are  negligable.  Usually  frequencies  less  than  lKHz  are  alright 
but  they  should  be  <  1Hz  to  ensure  no  induced  polarization  effects.  Alternatively,  one  can  use  exotic 
electrodes  such  as  Ag-AgCl  to  suppress  .  ontact  potential. 

Another  problem,  surprisingly,  is  the  high  sensitivity  of  the  methods.  Conductivities  span  a  wide  range 
and  these  directly  affect  the  measu*  potentials.  Thus  a  perfect  conductor  and  one  an  order  of  magnitude 
greater  than  the  “background"  conductivity  are  virtually  indistinguishable.  This  makes  it  very  difficult  to 
accurately  measure  conductivity  values  when  high  conductivity  contrasts  exist.  Electromagnetic  induction 
gets  around  this  since  the  eh':ct  of  including  the  magnetic  interaction  is  to  make  the  response  nonlinear 
in  <r. 


3.2  A  Note  On  Ohm’s  Law 

Fundamental  to  d.c.  resistivity  methods  is  Ohm's  Law: 

J  =  <r£  (121) 


where  J  is  the  current  density,  <r  is  the  scalar  conductivity  and  E  is  the  electric  field. 

Note  that  Ohm’s  law  is  linear.  This  is  fortunate  since  other'-’ «  most  problems  would  be  intractable. 
It  is  an  empirical  relationship,  however,  and  a  number  of  materials  do  not  obey  the  law  for  high  current 
densities.  For  most  experiments  this  is  not  a  problem,  since  |J|  <  lA/m1  except  near  the  electrodes. 
Still,  one  should  be  aware  of  the  limitations  of  the  approximation  and  approach  each  problem  cautiously. 

Equation  121  is  actually  the  scalar  form  of  Ohm’s  Law  and  is  true  only  for  &n  isotropic  medium.  In 
general,  W  is  a  rank-2  tensor  and 

f=WM  (122) 

To  show  that  anisotropy  can  be  significant,  corsider  three  simple  dual  conductivity  models  with 
conductivities  aj  >  <rj  (Fig.3.4): 

Model  (a)  ia  a  set  of  orthogonal  rods  of  conductivity,  cri,  imbedded  in  medium  of  conductivity,  <rj, 
In  (b)  the  rods  are  replaced  by  parallel  sheets  and  in  (c)  by  random  spheres.  For  all  model*  the  volume 

44  The  contrail!  are  nriiuuily  due  tf?  difference*  in  nil  ally  between  fluid*  internal  to  and  (unrounding  organ*. 
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fir  otion,  p,  of  material  1  is  much  less  than  2.  These  models  simulate  minerals  of  low  conductivity  with 
interstitial  water  or  senu-  metallic  inclusions.  The  effective  conductivities  are  found  to  be: 


cm  =  <r„  =  <rt  :=  +  (1  -  p)<r, 

(123) 

(b)  /T#  =  £Ty  =  p<ri  -f  (i  -  p)o>3 

<Ti  =  <Tl<T  if  (po-j  +  (1  -  p)(7l] 

(124) 

/'2a-2-rir1  +  2p[<ri-rra]\ 

:  CTy  =  0,  S3  — - - - ==• - r  <r3 

V  20-J  +  <Ti  -  p  £Tl  -  <TJ  / 

(125) 

To  illustrate  the  variation,  we  compare  the  results  of  the  three  models  in  Table  3.2.  It  is  assumed 
that  the  volume  fraction  p=  0.2  and  ff\  —  1 0<>-2  or  100<rj. 

We  ne'e  from  this  that  even  in  isotropic  media,  measured  conductivity  is  strongly  dependent  on 
geometry.  Also,  “platy*  structures  are  highly  aeolotropic  (anisotropic).  Aeolotropy  is  very  difficult  to 
analyse  and  generally  we  must  make  the  isotropic  assumption  if  a  mathematical  model  is  desired. 

3.S  Hock  Conductivity 

Assuming  that  electrostatic  measurements  can  be  made  and  the  inverse  problem  solved  to  ermine  the 
conductivities  of  the  body  being  measured,  one  must  attempt  to  determine  the  material’s  identity  from 
it*  conductivity.  When  imaging  the  human  body,  for  example,  this  may  not  be  too  difficult  since  bone 
has  quite  4  different  conductivity  than  muscle  or  fat.  Much  of  remote  sensing  is  related  to  geoexploration 
and  the  picture  for  rocfc/mincrai  conductivity  U  considerably  more  confusing. 

All  minerals,  except  the  smi-metallica,  are  insulators  with  10~lT  <  cr  <  10~11  S/pi.  Impurities,  such 
as  electron  donors/acceptors  and  other  crystal  defects,  increase  this  by  a  few  orders  of  magnitude.  In 
the  field,  however,  rock  conductivities  are  ~  1Q~S  to  10~l  S/m  (Table  3.1).  This  glarin'*  discrepancy 
is  accounted  for  by  the  fact  that  rock  conductivities  have  little  to  do  with  mineral  composition  but 
rather  are  related  to  permeability  /porosity  of  rock  to  fluids  and  the  conductivity  of  the  interstitial  fluid, 
Much  work  has  been  doni  and  many  models  have  bec»  developed  to  determine  suets  properties  fton 
measured  conductivities  but  it  is  beyond  the  scope  of  this  course.  There  are,  however,  general  trends  of 
conductivities  with  rock  type  and  these  are  presented  graphically  in  Fig.3.5.  Not  shown  in  f.g.S.fi  a*e  til* 
metallic  minerals.  These  include  free  metals  (<r  ~  10t+#-10+rS/m)  and  sullidc  ores  (a  ~  10+s~10+,*S/«r). 
Structural  conductors,  i.e.,  those  with  faults  cr  fracture  tones  that  nip* u re  large  quantities  of  water  are 
also  excluded, 

3*4  Theory 

3.4.1  General 

The  basic  idaa  of  resistive  measurements  is  to  employ  an  array  of  electrodes  connected  to  the  body  of 
interest  in  a  (hopefully)  appropriate  geometry  (Fig.S.fl).  The  quantity  of  interest  is  2  —  V/l  which  is 
called  the  “transfer  impedance”.  For  d.e.  or  near  d.c.  frequencies,  the  quantity  becomes  2  H,  the 
Transfer  resistance".  We  will  sec  that  for  most  geometries  that  can  be  simply  analysed, 

R  ~  Fp .  (120) 

where  p.  is  the  “apparent  resistivity"  of  the  body  and  F  is  a  geometry-dependant  quantity  which  has 
units  of  inverse  length.  We  wish  to  interpret  the  apparent  resistivity  in  terms  of  the  structure  of  the  body 
being  measured.  This  is  soldosa  straightforward. 

We  will  assume  in  the  following  discussions  that  the  medium  is  isotropic  and  thus  conductivity  is  a 
scalar  quantity.  Then  Ohm’s  law  is  given  by  Equation  121.  Furthermore,  we  assume  that  the  eiectrie 
field  is  conservative,  which  is  equivalent  to  assuming  that  we  are  in  the  very  low  frequency  regime.  Thus 

Vxlsd  (127) 

B  =  -Vj  (12$) 


24 


where  ip  is  the  scalar  electric  potential  Except  at  current  injection  sites,  the  time  rate  of  change  of  charge 
is  sero  and  hence  the  equation  of  continuity  becomes 

V  ■  /=  0  (129) 

Substituting  Equation  128  in  129  yields 

V  •  (<rVip)  =  0  (130) 

which  may  be  written 

VV+-(VV-)-(V<r)=:  0  (131) 

The  basic  Equation  130  or  131  together  with  the  appropriate  boundary  conditions  uniquely  specify 
the  problem  and  guarantee  a  unique  solution.  Assume  two  regions,  1  and  2,  and  an  interface  surface 
between  them  whose  unit  normal  at  a  point  is  n.  If  the  displacement  variable  along  the  normal  direction 
is  denoted  n,  then  the  boundary  conditions  are 

V>1  =  V1!  Dirichlet  b.c.  (132) 

Neumann  b.c.  (133) 

Note  however,  [5]  (pp.40-45),  that  a  solution  to  the  Poisson  Equation  130  with  ip  and  dip /On  specified 
arbitrarily  on  a  cloted  boundary  (the  so-called  Cauchy  boundary  conditions)  does  not  exist!  This  is 
because  only  a  single  boundary  condition  (Neumann  or  Dirichlet)  »a  necessary  at  each  point  of  the 
bounding  surface  to  define  a  unique  solution  to  the  Poisson  equation.  We  shall  use  this  as  a  method  of 
solving  the  inverse  problem  later. 

As  with  magnetostatics,  a  wide  variety  of  models  have  been  developed,  although  almost  paradoxically 
there  are  only  a  handful  of  general  geometries  which  yield  closed  form  solutions.  We  shall  solve  a  simple 
model  as  an  example  of  the  methods  to  use.  Although  the  example  is  not  a  compact  object  geometry,  it 
does  illustrate  typical  solution  methods.  The  object  of  the  problem,  once  again,  is  to  deduce  underlying 
anomaly  properties  by  resistance  measurements  on  the  surface  of  the  body.  The  first  step  is  to  have  the 
mathematical  model  and  the  second  is  to  solve  the  inverse  problem. 

A  case  of  importance  in  geophysics  is  one  in  which  <r  is  a  function  only  of  x.  Cylindrical  coordinates 
(r,  4>,  x)  are  then  a  good  choice  due  to  asimuthnl  symmetry  and  131  becomes 


d%ip  1 0*V>  1  8ipQ<r 

8r *  r  fir  fix’  ^  <r  fix  fix 

(134) 

This  equation  is  solved  by  the  standard  method  of  separation  of  variables. 

V»(r,x)  =  R(r)Z{x) 

(135) 

where 

d*R  1  dR  ...  n 
dr 2  r  dr 

(130) 

and 

d'Z  1  dedZ  , 

it1  <7  dx  dx 

(137) 

with  A  being,  so  far, 
written 

an  arbitrary  separation  constant.  A  general  solution  of  Equation  134  can 

fOO 

then  be 

V-(r,x)=/  F(X)R(X,r)Z(X,x)dX 

Jo 

(138) 

where  A  is  obviously  positive  definite  and  f  ( A)  is  chosen  to  fit  the  boundary  conditions. 
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3,4.2  Example  of  Electrostatic  Modelling 

Consider  now  a  two  layer  earth  model  in  which  the  current  is  injected  into  the  top  layer  at  the  origin 
of  the  coordinate  system.  The  conductivity  is  constant  in  each  layer  (Fig.3.7).  First  let’s  examine  the 
field  in  the  vicinity  of  the  current  electrode.  Assume  the  electrode  to  be  a  small  hemisphere  of  radius  b 
(Fig.3.8).8*  The  radial  current  density  is  given  bySfl 


J,= 


2 


where  a  total  current  I  is  injected  at  the  electrode  and 

Jis  =  2s  +  r* 


The  corresponding  radial  electric  field  is 
where 

near  the  electrode.  Thus  the  primary  potential  (due  to  the  current  electrode)  is 

i>P  = 

From  Bessel  function  theory  (i.e.  [10]), 

~  =  /  Jo(Ar)e-**dA  for  2  >  0 

H  Jo 


Er  =pi  Jr 


Pi  ~  <rx  1 


if-  I(h 


2t  R 


(139) 

(140) 

(141) 

(142) 

(143) 

(144) 


Jo  i®  the  Bessel  function  of  the  first  type  of  order  0.  Now,  solutions  of  Equation  136  are  cither  Jo  (Ar)  or 
Y0  (Ar),  the  latter  being  the  Bessel  function  of  the  second  type.  But  Yb  (Ar)  is  singular  at  r  =  0  and  must 
be  excluded.  Solutions  of  Equation  137  for  constant  <r  are  of  the  form  exp(±Ax).  Thus,  the  resultant 
potential  in  the  upper  layer  is 

i>i  =  i>p  + 

/oo 

[i(A)e-^  +  ^(A)e+^]  J0(Ar)dA  (145) 


for  0  <  2  <  h.  Equations  144  and  145  can  be  combined  to  yield 


* = r {[i + ^(x)]  e_** + b[x)**x'  }  Mxr)dx 


In  the  lower  layer,  since  z  -*  00,  the  potential  has  the  form 

=  ~  f"  c(\)*-XtM\T)d\ 

To  find  the  dimensionless  constants  A,  B,  C,  we  must  apply  boundary  conditions.  These  are 

-£-  =  0  at  2,0 

^1  =  V’J  nt  2  =  h 
Hi  dl(>3 

al  t=h 


(146) 


(147) 


(148) 


11  ll  is  assuratd  implicitly  that  the  second  electrode  U  si  infinity. 

**  Boundary  conditions  are  imposed  implicitly  la  this  equation  since  the  2r  (actor  ensures  that  all  current  must  flow  into 
the  s  >  0  region. 
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Applying  these  to  Equations  146  and  147,  we  obtain 


A-B-  0 


(1  +  A)  e~Xk  +  Be+Xk  =  Ce~xk 
<ri  [(1  +  A)  e~Xh  -  Pe+Afcj  =  <T2Ce~Xh 

(149) 

These  yield 

A~B  =  e~iXkK/  (l  -  Ke~3Xk) 

(150) 

where 

K  =  (pi  -pi)/(ps  +pi) 

(151) 

Thus,  the  solution  for  the  potential  in  the  upper  layer  is  Equation  146,  with  A, 
151.  Most  often,  measurements  of  potential  are  made  on  the  earth’s  surface  ( 

B  given  by  Equations  150, 
z  =  0).  Then 

(152) 

where 

f<x>  -J  xh 

<?(.,*)  =  l  +  2K,jf 

(153) 

The  limiting  cases  are: 

G(r,lif)— »1  as  hjr  — ►  oo 

G(r,K)  -*  —  as  h/r  -  0  (164) 


As  nn  interesting  (?)  aside,  by  expanding  the  denominator  of  Equation  153  in  a  power  series  in 
K  exp  (-2A h),  it  is  easy  to  show  that 


G(T,K)~l  +  2Kr'£iKnIn 


»=o 


where 


u=  re-«*+')»M\T)d\=- - i - -Hi 

Jo  (r*  4- <[« -f- 1]  2A}a) 


(165) 


(156) 


The  secondary  potential  is  thus  seen  os  being  due  to  an  infinite  series  of  image  sources  at  positions 
(0,0,  [n+  l]2h)  <md  strength  2JJC*+l. 

3.4.3  Four  Electrode  Array 

The  one  electrode  configuration  is  not  very  practical  and  so  we  consider  a  more  general  geometry  (Fig.3.0). 
Two  current  electrodes,  a  source  of  4/  amperes  and  a  sink  of -I  amperes  are  located  on  the  s  =  0  plane- 
The  voltage  difference  between  points  Pi  and  Pj,  also  on  the  z  =  0  plane,  Is  measured.  Without  loss  of 
generality,  we  can  place  the  source  electrode  at  the  origin  and  the  sink  on  the  x  axis.  By  superposition, 
the  potential  V*i  at  P.  is 


i>(* lift)  - 


Ipx 

2* 


*7(^111*0  _  Q(fsi)  K) 
ru  rji 


where 


'll  =  »i  +  V? 

rji  =  (*t  -  »o)*  +  V? 


The  potential  at  Pj  is  similarly 


y(Tl-n)-M[0(r..,K)  Q(r.„K)- 

2t  rn  fjj 


(167) 

(168) 

(159) 
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where 


ris  =  A  +  i/j 
r\i  =  (a2  -  a0)1  +  y\ 

The  voltage  between  the  electrodes,  V  =  fa  -  fa,  is 

G{rn,K)  G  (rji,  K)  G{rlitK)  {  G(r»,K ) 

^11  Tit  T32 


(160) 


(161) 


Now  for  homogeneous  ground  (pj  =  p\  or  h  — >  oo),  G  —  1  at  all  positions.  We  define  the  “apparent 
resistivity”,  p«,  as  that  resistivity  which,  for  a  homogeneous  earth,  would  yield  the  same  voltage  difference 
as  for  the  inhomogeneous  model.  Clearly, 


gfWO  _  O G(n„iC)  , 

_  Hi _ r,j  Tn  J  ral 

pi  “  -1 - 1 — —  -f  — 

r  »u  rai  ')>  'll 

For  the  Wenncr  array  (Fig.  3.2),  m  =  rn  =  a,  and  ni  =  r ai  =  2a,  and 


^  =  2  G(a,  K)-G(2a,K) 


(162) 


(163) 


This  function  is  shown  in  Fig.3.10.  We  see  that  p«  >  pi  for  p*  >  pi  and  p.  <  pi  for  pj  <  pi.  Also,  as 
the  upper  layer  thickness  becomes  much  greater  than  the  electrode  spacing,  p«  -*  pi ,  as  is  expected. 


3.5  The  Inverse  Problem 

Curves  of  the  type  just  shown  are  routinely  used  to  provide  a  rough  interpretation  of  resistivity  data 
in  electrical  depth  sounding  applications,  provided  that  the  underlying  structures  are  simple.  Vertical 
models  are  available  for  the  3-laycr  case,  2-layers  with  horisontal/ vertical  anisotropy  and  p  oc  exp  (-ax). 
Horisontal  models  are  available  for  the  vertical  dyke  and  vertical  wedge.  Other  models  are  available  for 
the  homogeneous  spherical  inclusion,  several  cylindrical  geometry  models,  and  spheroidal  inclusions. 

In  geophysics,  the  inverse  problem  is  usually  solved  merely  by  choosing  the  most  appropriate  geometric 
model  and  attempting  to  interpret  the  effective  resistivity  based  on  curves  such  as  Fig.3.10.  If  the  model 
is  wrong,  this  fails  badly  and  often  there  is  no  indication  that  it  failed.  Often,  there  is  no  appropriate 
model,  although  one  may  require  insight  from  other  remote  sensing  measurements  to  ascertain  that  the 
model  is  unsuitable. 

Because  of  the  limited  choice  of  models,  considerable  work  has  gone  into  solving  the  inverse  problem 
by  numerical  methods  which  can  handle  general  geometries.  Many  of  these  are  “impedance  tomography" 
techniques,  which  assume  that  the  current  in  the  body  of  interest  follows  ray-like  paths.  This  allows 
powerful,  well-established  techniques,  such  as  back- projection,  Radon  transform,  etc.,  which  are  used  in 
conventional  X-ray,  gamma,  PET  or  NMR  tomography  to  be  applied  to  the  problem.*7  There  are  a 
number  of  problems  with  methods  requiring  the  assumption  of  ray-like  current  paths.  First,  the  meth¬ 
ods  require  the  measurement  of  voltages  at  active  current  electrodes.  This  necessarily  introduces  the 
effects  of  contact  and  spreading  resistance  which  degrades  the  measurements.  This  is  a  major  concern, 
since  measurements  have  to  be  accurate  to  produce  fine  detail  in  the  tomographic  image.  Second,  the 
experimental  techniques  necessary  to  constrain  the  problem  to  the  ray-like  path  assumption  are  compli¬ 
cated,  usually  involving  complex  guarded  electrodes.  Third,  current  paths  are  highly  dependent  on  the 
object  of  interest,  in  particular  its  conductivity  distribution,  even  when  guard  electrodes  arc  used,  and 
often  the  ray  assumption  is  not  justified.  Since  it  it  the  conductivity  distribution  we  seek,  one  cannot 
know  a  priori  whether  the  ray  assumption  is  appropriate.  A  more  in-depth  discussion  of  the  problems 
associated  with  conventional  tomographic  techniques  it  found  in  [52].  In  spite  of  all  this,  a  number  of 
researchers  have  had  some  success  obtaining  impedance  tomography  images.  Henderson  and  Webster  [47] 
have  obtained  isoadmitlance  contours  of  the  chest  using  an  array  of  guarded  electrodes.**  Other  medical 

3TOne  doci  not  neceassxily  have  to  assume  straight-line  paths  for  some  of  the  more  sophisticated  tomography  algorithms. 
For  example,  some  researchers  have  assumed  curved  current  flux  tubes. 

41  The  guarding  was  intended  to  force  a  long  narrow  measurement  volume  centered  on  each  electrode. 
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applications  and  methods  are  discussed  in  [51],  geophysical  applications  may  be  found  in  [48]  and  ground 
water  pollution  monitoring  is  discussed  in  [56]. 

An  alternative  approach  is  to  inject  current  at  one  or  more  sites,  measure  the  surface  potential  distri¬ 
bution  and  numerically  solve  the  field  equations.  The  inverse  problem  associated  with  a  single  excitation*8 
does  not  generally  have  a  unique  solution.  By  this,  we  mean  that  different  internal  conductivity  distribu¬ 
tions  can  give  rise  to  the  same  surface  potential.40  If  a  number  of  excitations  are  used  and  the  potential 
distribution  measured  for  each  excitation,  the  indeterminancy  can  be  reduced.  Each  excitation  is  akin  to 
a  different  view  in  conventional  tomography  and  hence  it  might  seem  intuitive  that  if  a  sufficient  number 
of  independent  excitations  are  used,  the  conductivity  distribution  might  be  estimated  with  reasonable 
certainty.  Since  knowledge  of  the  conductivity  distribution  is  effectively  an  image  of  the  conductivity ,  we 
will  call  the  method  “conductivity  imaging”.41  The  development  closely  follows  that  given  by  Wexier  et 
al.  ([54],  [55]). 

The  geometry  is  shown  in  Fig.3.11.  An  array  of  electrodes  is  placed  on  the  surface  of  the  body  to 
be  measured.  Current  is  injected  at  one  electrode  and  extracted  at  another.  Potential  measurements 
are  made  at  the  other  electrodes.  These  are  referenced  to  the  potential,  <f> rbF)  at  an 
arbitrary  but  fixed  electrode.  The  latter  process  constitutes  one  excitation.  Note  that  potential  measure¬ 
ments  are  not  made  at  current  sites.  The  electrode  array  shown  is  merely  illustrative.  The  actual  one  used 
may  not  be  square  and  may  consist  of  more  or  less  electrodes,  depending  upon  the  actual  experiment. 

The  resulting  surface  voltage  distribution  is  measured  for  several  sets  of  excitations  and  stored  for  later 
processing.  The  measurements  may  be  taken  using  an  automated  digitally  controlled  data  acquisition 
system.  The  general  method  we  will  employ  is  to  guess  at  the  conductivity  distribution  and  use  it  to 
calculate  a  potential  distribution  throughout  the  volume  and  in  particular  at  the  surface.4*  The  guess  is 
unlikely  to  agree  with  the  true  subsurface  conductivity  distribution  and  so  the  calculated  and  measured 
surface  voltages  will  disagree.  An  algorithm  which  will  iteratively  refine  the  conductivity  is  needed. 
This  will  then  be  repeatedly  applied  until  the  measured  and  calculated  surface  voltages  are  in  what  is 
considered  to  be  acceptable  agreement. 

The  equations  governing  the  problem  may  be  easily  derived.  If  we  assume  that  the  local  conductivity 
is  denoted  k  and  the  frequency  is  low  enough  so  that  the  electric  field  E  may  be  expressed  in  terms  of  a 
scalar  electric  potential  i>, 

i  -  -Vi,  (104) 

Then  by  the  equation  of  continuity 

V./=-|  =  -/  (165) 

where  J  is  the  current  density  and  /  is  the  impressed  current  source  distribution  within  the  volume  of 
interest.  These  two  equations  with  Ohm’s  Low  yields  the  Poisson  equation  for  continuous  inhomogeneous 
media4* 

V’*Vi>  =  -f  (106) 

For  ease  of  computation,  the  region  of  interest  (Fig.3.11)  is  assumed  to  be  bounded  by  the  measurement 
surface  and  five  orthogonal  plane  faces.  If  the  measurement  surface  were  a  horisontal  plane,  these  six  faces 
would  form  a  cube.  The  volume  of  interest  is  divided  into  a  finite  element  grid.  In  the  simplest  scheme, 
potentials  are  computed  for  each  excitation  at  the  node  points  (mesh  intersections)  and  the  conductivity 
iB  then  estimated  within  the  intervening  regions.  The  conductivity  distribution  in  the  region  of  interest 
is  initially  assumed  to  be  uniform  and  the  iterative  procedure  to  improve  the  conductivity  estimates  is 
then  applied.  The  procedure  is  outlined  below: 

11  The  injection  of  a  fixed  current  si  a  specific  site  and  withdrawal  from  another  ipeciflc  titc  conilltutea  an  excitation. 

*°Thii  war  acen  in  the  2-layer  earth  model  (iee  Fig.3.10),  Alio,  consider  the  simple  example  of  an  array  of  electrode* 
on  the  earth  •  surface,  all  at  the  tame  potential.  Thii  implies  that  a  conductor  it  immediately  beneath  the  measurement 
turf  ace  but  do*s  not  lay  whether  the  conductor  ia  a  thin  sheet  or  extends  to  large  depths. 

41  This  can  lead  to  some  confusion  since  the  method  ia  also  referred  to  in  the  literature  as  impedance  tomography.  In  this 
section  "impedance  tomography"  will  be  reserved  for  those  techniques  assuming  ray-like  current*. 

41 A  numerical  method  such  as  the  finite  element  method  is  employed. 

44  Equation  1M  is  strictly  true  only  for  aero  frequency.  For  slightly  higher  frequencies,  it  ia  valid  if  the  conductivity  is 
treated  as  a  complex  quantity  (sec  next  section).  For  still  higher  frequencies,  the  Helmholta  equation  must  be  used. 
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1.  Calculation  of  potential  and  current  density  with  Neumann 

boundary  conditions  -  Given  the  generally  inhomogeneous  conductivity  distribution  from  the 
previous  iteration,  the  potential  distribution  is  computed  by  solving  the  Poisson  Equation  166. 
The  potential  distribution  is  in  turn  used  to  calculate  the  current  density  distribution  J  (A/m3). 
This  latter  quantity  is  what  we  actually  seek. 

In  solving  for  the  potential,  the  inhomogeneous  Neumann  boundary  conditions  ore  used,  namely 

«C)(tcL“*W 

where  t  denotes  a  point  on  the  surface  bounding  the  volume  of  interest,  subscript  “1”  denotes  a 
quantity  derived  from  the  first  step  of  the  iteration,  n  is  the  displacement  in  the  direction  normal  to 
the  surface  and  h(»)  (A/m3)  is  the  current  density  entering  or  leaving  the  medium  over  the  element 
of  bounding  surface  at  that  point.  The  quantities  h(»)  are  derived  from  the  measured  currents 
at  the  electrodes.44  Where  no  current  is  impressed,  h(t)  =  0.  This  is  true  for  all  points  on  the 
bounding  surface  other  than  the  current  input  and  output  electrodes.  The  justification  for  this  on 
the  measurement  surface  is  that  the  conductivity  of  air  is  approximately  0.  For  the  other  faces,  it 
is  implicitly  assumed  that  the  bounded  volume  is  sufficiently  large  that  negligable  current  crosses 
those  faces. 

The  boundary  conditions,  assumed  conductivity  distribution  and  Poisson  equation  are  then  used 
to  calculate  the  potential  distribution  throughout  the  volume,  using  one  point  on  the  measurement 
surface  as  a  reference.  This  is  normally  done  by  a  finite  element  method  which  computes  the  field 
by  dividing  the  region  into  numerically  manageable  intervals  called  “finite  elements".  The  method 
calculates  the  potential  at  the  node  points,  using  one  of  a  number  of  field  calculation  methods 
(i.e.,  [53]).  In  simple  versions,  the  nodes  are  vertices  of  the  cube  elements,  but  more  sophisticated 
versions  allow  for  nodes  to  be  at  other  points.  After  the  potentials  are  found  at  the  node  points, 
located  at  mesh  intersections,  the  potentials  at  other  points  may  be  obtained  by  interpolation. 

Next  we  find  the  electric  field  throughout  the  volume  from  Equation  164  and  then  the  current 
density  distribution  from 

Jl  =  kEi  (168) 

Step  1  is  repeated  for  all  other  excitations  and  the  ipsultii.3  fields  are  stored.  41 

2.  Calculation  of  potential  with  Dirichlet  boundary  conditions  -  Next  the  interior  potentials 
are  calculated  again  by  an  independent  method,  namely  by  the  use  of  the  Dirichlet  boundary 
conditions, 

M»)  =  S{*)  (160) 

on  the  measurement  surface.  (Here  the  subscript  “2"  denotes  quantities  derived  from  the  second 
step  of  the  iteration,  and  g($)  are  the  potentials  measured  at  the  bounding  surface,  namely  the 
measured  fa  at  the  potential  electrodes.)  In  addition,  the  boundary  conditions  must  include  the 
Neumann  condition  at  the  current  electrodes  since  measured  potentials  are  not  available  at  these 
sites.4®  Boundary  condition  at  the  reference  electrode  is,  of  course,  g(t)  ~  0  and  the  Neumann 
condition  h(t)  =  0  is  necessary  on  the  other  five  faces.  The  quantity  we  wish  to  derive  from  this 
step  is  V*s- 

3.  Calculation  of  conductivity  -  As  was  stated  in  Section  3.4.1,  the  Poisson  equation  yields  a  unique 
solution  \!>  for  specified  k,  when  a  single  boundary  condition  (Neumann  or  Dirichlet)  is  specified  at 
each  boundary  point.  In  fact,  the  solution  V>  can  be  used  to  derive  the  other  boundary  condition  at 

44  The  integral  of  h[i)  over  the  electrode  surface  yields  the  total  current  /.  Electrodes  can  be  modelled  as  hemispheres  or 
points,  the  latter  often  being  convenient. 

4>The  current  distribution  will  initially  be  approximate  because  the  conductivity  and  hence  the  potentials  are.  However, 
reasonable  current  Cow-line  patterns  an  obtained  even  for  very  approximate  ft.  This  Is  because  the  current  is  constrained 
to  enter  the  measurement  surface  through  one  electrode,  then  spreads  widely  through  the  volume  and  must  converge  toward 
the  other  current  electrode. 

44  Voltages  measured  at  the  current  Injection  electrodes  an  prone  to  error  due  to  contact-potential  problems  (see  Section 
3.1). 
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eaeh  point.  What  this  means  then,  is  that  if  the  actual  conductivity  distribution  k  is  known,  Ohm’s 
Law  states  that  Ji  =  nVrj) j.  But  the  boundary  conditions  are  based  on  measurements,  subject  to 
error,  and  k  is  an  estimate  and  thus  the  two  different  sets  of  boundary  conditions  will  not  yield  the 
same  internal  fields. 


We  can,  however,  adjust  k  to  ensure  that  the  two  sets  of  boundary  conditions  are  as  compatible 
as  possible  in  some  sense,  by  minimizing  a  quantity  related  to  some  average  over  all  points  and 


excitations  of 
minimize 


fi  +  KV^ijj.  If  we  seek  to  optimize  the  solution  in  a  least  squares  sense,  we  should 

R  =  X)  /  /  Jv  (£  +  rt'h)  •  (ji  +  KVfa)  dv  (170) 


where  R  is  the  squared  residual  sum,  V  is  the  region  of  interest  and  X  denotes  the  excitations  used. 
Because  we  have  employed  a  finite  element  scheme,  the  integral  over  V  is  really  a  sum  of  integrals 
over  the  element  volumes  Vj,  i.e., 


^  =  E  E  /  /  jy{fl+  K^3)  ‘  +  K^3)  ^ 

=  ??///,  ^1  •  Ji  +  2kj  Ji  •  VV>a  +  dv  (If  1) 

where  Kj  is  the  conductivity  distribution  within  element  j. 

The  size  of  the  elements  may  be  chosen  so  that  kj  is  constant  within  each  element  and  then  we 
minimize  the  residual  by  demanding 

~  =  E  /  J  Jy  (2*  • **  +  ■  V^a)  dv  =  0  (172) 


where  and  V>a  are  fixed  at  the  previously  estimated  values.  We  can  rearrange  the  last  equation 
to  give  an  estimate  of  the  optimum  conductivity  distribution  for  the  iteration, 


E.///v4  VV>a- V^adv 


(173) 


This  last  equation  is  then  applied  over  all  elements  in  the  volume  to  obtain  the  revised  estimate  of 
the  conductivity  distribution.47 


4.  Recursive  improvement  -  Steps  1,  2  and  3  constitute  an  iteration.  The  next  iteration  begins 
again  by  solving  the  Neumann  boundary  value  problem  for  all  excitations  using  the  new  conduc¬ 
tivity  distribution  estimate.  The  computed  potentials  at  the  boundary  are  compared  with  the 
measured  ones.  If  the  differences  exceed  some  a  priori  thresholds  or  if  insufficient  iterations  have 
been  performed4*,  the  iteration  continues  with  the  solution  of  the  Dirichlet  boundary  value  prob¬ 
lem.  Otherwise  the  resulting  image  (2D  or  3D  depending  on  the  problem)  of  the  conductivity 
distribution  is  processed  by  one  of  a  number  of  standard  techniques  and  is  output. 

As  a  simple  but  somewhat  contrived,  example  of  the  three  dimensional  case,  a  cube  with  four  layers 
of  finite  elements  is  shown  in  Fig.3.12  ([54]).  Layers  two  and  three  contain  a  square  object  whose 
conductivity  is  five  times  the  host  medium.  Simulated  “measurements”  have  been  made  at  the 
top  surface  only.  Clearly  the  estimation  improves  with  iteration  count.  Initially  a  conductivity 
artifact  appears  at  the  surface  but  it  later  disappears.  Note  that  a  large  number  of  iterations  are 
needed  and  even  then  the  image  is  quite  crude.  Objects  of  a  more  general  shape  also  tend  to  be 
crudely  imaged.  Furthermore,  computer  time  is  excessive,  being  of  the  order  of  an  hour  on  a  large 
mainframe  computer.  Work  is  underway  to  solve  both  the  accuracy  and  speed  problems. 

4TNote  that  no  malrir  operation!  are  Decenary. 

41  The  minimum  necewarjr  number  of  iteration!  U  a  heuristic  parameter  which  la  derived  after  considerable  user  experience 
on  a  particular  problem. 
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3.6  Induced  Polarization  (IP) 

Assume  now  that  rather  than  d.c.,  we  are  operating  at  a  very  low  frequency,  u.  Again,  as  for  electrostatics, 
we  assume  that  magnetic  effects  can  be  neglected.  Frequently  we  find  that  <r  is  a  function  of  u>  and  is  a 
complex  quantity,  i.e.  the  current  is  out  of  phase  with  the  voltage.  Thus,49 

a1  (iu)  =  a  (iu)  +  iue  (iw)  (174) 

with  <r,  e  real  and 

f  (iu)  =  a>  (iu)  E  (iu)  (175) 

Components  of  J  and  E  are  phasors  and  thus,  for  example,  the  real  physical  quantity  e«(t)  (subscript  z 
denotes  the  z  component)  is  related  to  the  complex  phasor  by 

e.(t)  =  St  {^(iuOe'"4}  =  l^.tcos  (ui  +  fa.)  (170) 

and 

\J.\  =  \*\\E.\ 

(177) 

where  <pj.  is  the  phase  of  J,,  <pa  is  the  phase  of  a  and  <)>e.  is  the  phase  of  Em. 

Problems  in  IP  are  usually  solved  by  assuming  a  model  for  p  =  (er1)-1  and  then  solving  a  particular 
boundary  value  problem  as  for  d.c.  conductivity  methods.50 
A  simple  model  for  p  is  the  JV-pole  model 


<’<•)  =  X>;~j - 

n=0  n 

with  an,  An  real  and  s  being  the  Laplace  variable.  Another  model  commonly  found  in  the  literature  is 

*)-Sau£st*- 

with  Gn,  Cn,  pa,  all  real.  This  allows  |p|  to  be  nonsero  as  «  -♦  oo,  The  Cole-Cole  model  is  also  very 
popular: 

P(*)  =  Po  (l  -  mo  ^1  —  [l +  (iwT)*]-1)}  (180) 

where  po,  mo,  r  are  real  constants  and  if  is  an  empirically  adjusted  noninteger. 

All  these  are  phenomenological  models  which  in  some  manner  describe  the  electrical  properties  of  the 
material  being  probed.  The  cause  of  IP  is  not  known  in  detail. 

Intuitively,  however,  if  the  impedance  is  complex  and  decreases  with  frequency,  it  must  be  related 
to  capacitive  effects  in  the  medium  (recall  that  displacement  current,  although  small,  has  not  been 
neglected).  These  capacitive  effects  are  due  to  electrochemical  and  other  electromagnetic  properties  of 
the  medium.  To  illustrate  this  we  use  a  simple  model: 

Assume  a  single  sphere,  resistivity  pi,  radius  a,  in  a  uniform  medium  of  resistivity  p  and  assume  there 
exists  a  (complex)  interface  impedance  tjm  (ohra-m*).  Primary  field  f?o  is  uniform  and  along  the  polar 
axis  of  a  spherical  coordinate  system  (r,  0,  <j> )  centered  on  the  sphere  (Fig. 3. 13).  It  is  assumed  that  u  is 
small  and  Laplace's  equation  holds.  The  analysis  is  routine.  Suitable  potentials  are: 

i>i=Aorcoa$  for  0<r<o  (181) 

=  -Eqt  cos  6  +  Ar~i  cos  0  for  r>o  (182) 

*'Thi*  it  really  the  tame  at  lumping  the  small  effect  of  dl (placement  current  into  an  “effective  conductivity". 

*°  Thlt  it  possible  since  the  frequency  is  assumed  sufficiently  low  so  that  the  Poitaon/Laplace  equations  apply. 


(179) 
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where  Ao,  A  are  constants.  Boundary  conditions  at  r  =  a  are: 

1  dij>  _  1  dij>i 
p  dr  pi  dr 


Vm  dlj) 


Equation  183  derives  from  J,  being  continuous  while  184  is  due  to  the  voltage  drop  across  r  —  a  being 
i}mJr.  Solving  Equations  181  to  184  yields 

j>  —  -EorcosO  +  a’^Eor-2  cos#  (185) 


*1  +  2* 
and 

*  =  —  +  —  (187) 

p  pa 

Now  assume  N  such  particles  (polar  axes  aligned)  in  a  spherical  region  of  radius  r0  (Fig.3.14)  and  if 
interactions  between  particles  are  ignored  then  at  point  P  we  have 


i  n  n  V —  EoaSX  ,, 

rfr  =  -Eor  cos  0  +  >  — cos 

i  r< 

Typically,  ro  <  r,  «  r  and  Oi  ss  6.  Then 

rl>  w  —Eor  cos 6  +  NaPxEor'*  cos  6  (188) 

But  if  we  assume  the  spherical  region  is  a  continuum  with  effective  resistivity  p„  then 

»  -25o*,cos0  +  —EocqsO  (189) 

*  P  +  *P* 

Equating  188  and  189  yields  the  important  result:51 

L  =  ±^2.  (i9o) 

P.  1  +  21.*  1  ' 

where  v  =  Nas/ro  is  the  volume  fraction  of  particles. 

If  we  assume  p\jp  <£.  1  and  v  <  1  (i.c.,  many  small  highly  conducting  particles  in  an  electrolyte)  and 
if 

6- — *  (191) 

2(iwr)*  V  ' 

then  we  obtain  the  Cole-Cole  equation  180. 

Equation  191  implies  that  the  interface  impedance  is  capacitive  in  nature.  The  Cole-Cole  model  has 
been  widely  used  since  it  fits  a  large  body  of  accumulated  data  quite  well. 

Finally,  it  should  be  noted  that  IP  measurements  are  commonly  made  in  the  time  domain.  In  such 
cases  it  is  assumed  that  the  time  dependent  impressed  current,  js(t),  is  known  (we  restrict  ourselves  to 
the  z  direction)  and  we  wish  to  predict  e,(t)  given  a  resistivity  model.  As  an  example,  if  we  assume  a 
model  for  p(t)  as  in  Equation  178M  and  if  we  assume  a  step  function  impressed  current, 

i.(t)  =  WO  (182) 


M •)  =  Jo/» 

*lThl»  U  analogous  to  the  Clautlua-MotoUt  equation  in  dielectric  theory. 
M  We  un  capital  letter*  to  denote  Laplace  variable*. 
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e. 


(194) 


"  r  i  i 

»(f)  =  Ii  1  (p(s)Ja(*)]  —  Jo  ^2  A-nL  1  /  ,  \ 

n=0  L'l'  +  OnJJ 

which  becomes 

e.(t)  =  JroE“[1-c"a“‘Mt)  (195) 

n=o  a* 

Letting  eo«  =  e,(oo),  we  can  plot  Am(t)  =  eB(t)/e o«  (normalised  step  frequency  response)  versus  i.  This 
is  shown  in  Fig.3.15. 

In  practice,  it  is  often  convenient  to  use  a  continuous  on-off  current  waveform  such  as  Fig. 3. 16.  The 
analysis  is  similar  to  above.  The  waveform  is  a  sum  of  shifted  step  functions  and  one  usually  analyzes 
the  “steady-state  transient  response”.  The  output  then  looks  like  that  of  Fig-3.17. 

Time  domain  induced  polarization  has  an  advantage  over  the  frequency  domain  in  that  many  fre¬ 
quencies  are  excited  by  one  measurement,  in  principle  providing  more  information.  Unfortunately,  the 
analysis  is  much  more  complicated  for  the  time  domain. 

4  Electromagnetic  Induction 

4.1  Introduction 

The  method  of  electromagnetic  induction  is  widely  used  for  detecting  conductive  objects  in  geophysir- 
nondestructive  testing,  detection  of  mines  and  artillery  shells,  archeological  exploration  and  treasure 
hunting.  A  conductive  object  is  exposed  to  a  time-varying  magnetic  field  (usually  in  the  Hs  to  KHz 
range)  produced  by  a  current-carrying  primary  coil.  The  resulting  eddy  currents  induced  in  the  object 
produce  a  secondary  time- varying  magnetic  field  which  is  detected  by  another  for  occasionally  the  same) 
coil. 

The  receive  coil  can,  of  course,  detect  variations  in  the  geomagnetic  field  and  it  is  important  that 
the  signal  due  to  the  geomagnetic  field  be  substantially  less  than  that  due  to  the  object  of  interest.  The 
geomagnetic  spectrum  is  shown  in  Fig.4.1.  Calculation  of  the  geomagnetic  “noise'*  is  dependent  upon  the 
receiver  design,  but  we  can  obtain  a  rough  estimate.  If  we  assume  a  receive  coil  of  radius  r,  a  geomagnetic 
flux  density  at  u>  Hs  of  2?”,  then  the  induced  emf  |  KJ1 1  at  frequency  w  is  given  by 

ICI  =  |*r*«S:3|r|  (196) 

where  the  bar  denotes  an  average  value  and 

S2  =  (197) 

Assuming  r  =  0.2m,  10“eT,M  u  =  2rx  10*Hs,  SJ  =  1/30  (from  Fig.4.1),  we  obtain 

jVjoool  w  2.6  x  10-8  Volts  at  1  KHs 

Sensitive  instrumentation  can  and  routinely  does  measure  such  small  voltages  in  the  course  of  fundamental 
studies  of  the  geomagnetic  field.  Note  that  noise  voltage  is  proportional  to  the  receiver  coil  cross  sectional 
area.  Geometries  associated  with  mines  and  artillery  shells  typically  use  coils  with  radii  of  less  than  1 
meter  and  detected  signals  are  of  the  order  of  1  x  10“aVolls  or  greater  so  geomagnetic  noise  does  not 
pose  a  problem.  In  nondestructive  testing  the  signal  of  interest  is  generally  much  larger  still,  owing  to 
the  small  coils  (radii  arc  in  the  mm  to  cm  range)  and  close-coupled  geometry  of  coils  and  object,  and 
geomagnetic  noise  is  even  less  of  a  problem.  In  geoexploration,  on  the  other  hand,  received  signals  are 
usually  very  weak  which  necessitates  large  receive  coils.  Geomagnetic  noise  must  then  be  contended  with 
and  occasionally  can  even  prohibit  exploration.  In  any  case,  coherent  averaging  can  increase  the  S/N  by 
substantial  amounts  since  the  geomagnetic  signal  is  essentially  stochastic. 

The  origin  of  the  geomagnetic  spectrum  is  beyond  the  scope  of  this  course  and  the  student  is  referred 
to  (60]  for  further  details.  We  only  note  that  the  effect  of  geomagnetic  noise  generally  tends  to  decrease 
with  increasing  frequency.  This  is  because  |V*j  is  proportional  to  w  and  the  spectral  density  decreases 
faster  than  w-1  with  the  exception  of  the  gyromagnetic  resonance  region. 

To  be  ttriclly  correct,  the  meaturemcnl  band  muil  be  ipedfled  in  order  to  determine  the  noire  at  a  particular  frequency. 
However,  thii  value  It  roughly  correct  for  typical  magnetometer  bandwidth!. 
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4.2  Principles  of  Detection 

4.2.1  General  Principles 

There  are  basically  only  two  types  of  electromagnetic  induction  detectors  -  continuous  wave  (CW)  and 
transient  detectors. 

Continuous  wave  detectors  generally  employ  a  single  frequency  sinusoidal  transmitter  waveform.  They 
may  then  be  further  subdivided  into  frequency  shift  detectors  and  mutual  inductance  detectors.  Frequency 
shift  detectors  usually  employ  a  resonant  oscillator  circuit  with  a  coil  acting  as  the  inductive  element. 
In  isolation,  the  oscillator  operates  at  a  fixed  frequency  determined  in  part  by  the  inductor.  If  the  coil 
is  now  brought  near  a  conducting  or  magnetic  body,  the  effective  inductance  changes  and  hence  the 
frequency  changes.  A  variety  of  techniques  are  employed  to  detect  the  frequency  shift.  One  of  the  most 
straightforward  is  to  mix  a  local  oscillator  signal  at  the  fundamental  frequency  with  the  shifted  frequency 
and  then  measure  the  beat  frequency  with  a  counter. 

Mutual  inductance  detectors  rely  instead  on  two  or  more  coils  in  close  proximity.  The  coils  are  oriented 
such  that  ideally  the  mutual  inductance  between  them  is  sero.  One  coil  (the  transmitter  or  primary)  is 
excited  by  a  single  frequency  sinusoidal  signal.  In  the  absence  of  an  object,  the  signal  coupled  into  the 
other  coil(s)  (receiver  or  secondary)  is  sero.With  a  conductive  or  magnetic  object  present,  the  mutual 
coupling  between  the  coils  changes  and  a  signal  is  induced  in  the  receiver  coil(s).M 

Transient  detectors  also  use  a  transmitter  coil  and  one  or  more  receive  coils.  Cut  rent  in  the  transmitter 
coil  is  switched  off  and  on  rapidly.  Transient  currents  are  induced  in  the  receive  coi!  (s)  even  in  the  absence 
of  a  target  object.  To  eliminate  this  mutual  coupling  signal,  some  systems  use  two  coaxial  receive  coils 
which  are  equidistant  from  the  transmitter  coil.  The  coils  are  wound  in  opposition  and  fed  to  a  difference 
amplifier  which  accordingly  produces  a  null  output.  A  target  placed  near  one  receive  coil  will  induce  a 
bigger  signal  in  that  coil  than  the  other  and  hence  a  signal  will  appear  at  the  amplifier  output.  The 
problem  with  this  method  is  that  the  signal  induced  by  the  primary  coil  in  the  secondary  coils  is  up  to 
10e  times  that  induced  by  the  object.  Thus  extremely  accurate  coil  matching  and  alignment  (balancing) 
is  necessary  if  the  object  is  to  be  detected. 

An  alternative  transient  approach  is  to  use  a  pulsed  transmit  current  and  to  only  measure  the  volt¬ 
age  in  the  receive  coil  during  the  quiescent  period  between  pulses.  Provided  the  transients  due  to  pri¬ 
mary  /secondary  mutual  coupling  decay  faster  than  those  due  to  the  target,  the  method  will  work  without 
the  need  for  critical  coil  balancing.  This  method  it  often  referred  to  as  “pulse  induction".  The  chief  draw¬ 
back  of  pulse  induction  over  the  balanced  3  coil  arrangement  is  that  pulse  induction  must  have  a  duty 
cycle  less  than  100  %  whereas  this  is  not  a  constraint  for  the  latter  method.  The  lack  of  a  need  for  critical 
coil  balancing  for  pulse  induction  more  than  makes  up  for  this  disadvantage. 

As  one  might  expect,  the  continuous  wave  and  transient  methods  are  connected  via  the  Laplace 
Transform.  We  shall  discuss  this  point  later.  Continuous  wave  detectors  generally  require  less  power  for 
a  given  sensitivity  than  do  transient  detectors.  On  the  other  hand,  CW  detectors  yield  less  information 
about  the  object  since  they  operate  at  a  single  frequency,  rather  than  exciting  with  a  broad  spectrum  as 
do  transient  detectors.  Furthermore,  mutual  inductance-type  detectors  require  very  critical  coil  alignment 
too,  which  may  detract  from  the  robustness. 

Many  detectors,  as  mentioned  previously,  use  a  transmit  coil  and  one  or  more  receive  coils.  The  choice 
of  geometries  is  almost  limitless,  although  they  fall  into  a  few  general  categories.  The  basic  idea  is  to 
have  a  mutual  inductance  that  is  sero  in  the  absence  of  a  target  and  thai  is  measurable  with  the  object 
present.  This  is  not  always  very  simple  to  do.  To  understand  better  how  to  achieve  this  goal,  we  must 
first  digress  to  discuss  the  field  of  a  coil  and  mutual  inductance  in  more  detail. 

4.2.2  Field  of  a  Coil 

Assume  a  loop  carrying  current  I  situated  in  the  x  -  y  plane  (Pig.4.3).  We  wish  to  calculate  the  vector 
potential  at  a  point  P,  which  because  of  axial  symmetry  may  be  chosen  without  loas  of  generality  at 

**  Other  variant*  um  a  circuit  which  very  new  owaUatlott.  The  tstkll  change  In  mutual  Inductance  induce*  oactUalloo. 
This  ha*  the  advantage  of  low  power  coaawnptloa. 
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<f>  =  0.55  The  vector  potential  is  given  by:54 


A(z) 


-  i  /  ZGU 


4  *Jv 


v  »<\ 


The  current  density  has  only  a  <t>  component  and  is  given  by 

J  =  —  J*sind>'z  4-  cos  <p'y  —  J$4> 

srwith 

J+  —  1  sin  (cos  6')  - — 


(198) 


(199) 


(200) 


Since  we  are  observing  at  <p  —  0,  the  2  component  of  the  current  density  will  yield  no  contribution  on 
integration.  Also,  since  J  has  only  a  4>  component,  so  does  A,  i.e,  A  =  A$4>.  Thus 
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Now  1/ 
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J*  "  *  , 

can  be  expanded  as  a  sum  of  spherical  harmonics,  K(m  (0,  <p). 
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(202) 


where  r<  (r>)  is  the  smaller  (larger)  of  a  and  r.  Subatitution  of  Equation  202  into  201  and  application 
of  the  properties  of  spherical  harmonics  yields 


s*wherc  Pj™  (cos  0)  is  the  associated  Legendre  function  and 


(203) 


(2n  -  1)!!  =  (2n  -  1) .  (2n  -  3)  •  -5  •  3  •  1 
Using  B  =  V  x  A  (sec  Section  4.3)  we  find  that 

A  =  t^  L_LL._L 

where  Pn  (co «0)  is  a  Legendre  Polynomial  and 

_  p/o>  ^  (-1)"  (2n  *  1)»  /2n+  2\  1  (r^  , 

~  ~T  ^  2*TnTi)!~  UTTJ  ?  U  {c0t6) 
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p/o^t-l)*  (2n-fl)ll  1  /ox*-  , 

4  ^  2-  (n  4-  1)!  r*\W  ,n+* 
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(206) 


where  the  first  equation  for  li$  is  for  r  <  a  and  the  second  u  for  r  >  a.  There  arc  two  important  cases. 
For  r  >  o  only  the  a  =  0  term  is  significant. 


u  _  co*  ® 

tfr  —  — ~ 

2  r3 

MAg«in,  primed  coordinate*  ar*  mure*  coordinate*,  u*pri«*d  ut  field  coordinate*. 
**TgU  will  fee  derived  in  Set  lion  S.3. 

,TWhy  hate  mt  feotkeeed  to  fcmIk  J  into  t  ud  p  cempMutli? 

M  Can  yon  dote*  l St  step*  tr oa  Equation  203  to  203? 
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These  are  just  the  field  components  of  a  dipole.  Thus,  sufficiently  far  from  a  loop  of  wire,  the  field  is 
dipolar  in  nature. 

For  r  o, 

..  T 

(208) 
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In  other  words,  close  to  the  coil  we  have  a  uniform  field  in  the  direction  of  the  coil  axis. 

4.2.3  Mutual  Inductance  of  Two  Coils 

We  are  now  in  a  position  to  discuss  the  mutual  inductance  of  two  coils.  We  assume  the  transmitter  and 
receiver  are  two  simple  loops  of  negligable  cross-section  and  are  situated  in  free  space.  If  a  sinusoidal 
current  flowing  through  loop  1  (the  transmitter)  is  denoted  Ij,  then  the  induced  emf  in  loop  2,  Vj,  is 


Vi  =  —iuMnh 


(209) 


where  u  is  the  angular  frequency  of  excitation  and  Mu  is  the  “coefficient  of  mutu.d  induction”  or  simply 
the  “mutual  inductance”.  By  the  reciprocity  theorem,  we  can  interchange  transndtter  and  receiver  and 
hence 

Mu  =  Mn  (210) 

Referring  to  Fig.  4.4,  the  magnetic  potential  at  point  aij  of  the  receiver  loop  is  given  by  Equation  198 
in  the  form 


_  po I  f  dlx 
4*'  fx  1*2  -  »1 1 


(211) 


where  subscript  “1”  refers  to  quantities  measured  at  the  transmitter  and  subscript  “2”  refers  to  quantities 
measured  at  the  receiver.  The  magnetic  flux  intercepted  by  the  receiver  is 


But 

which  implies  that 
Thus  we  have, 


$a  =  J  B  (»j)  •  d$i 
iVxI(aTj)]  -dst 

Vi  =  -tw$j 

Mu  =  $j/Ji 
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(214) 

(216) 


This  general  equation  can  be  very  difficult  to  solve  for  arbitrary  geometries.  For  the  special  case  when 
the  interloop  distance  greatly  exceeds  the  coil  diameters,  the  field  of  the  transmitter  at  the  receiver  is 
that  of  a  dipole.  The  general  expression  for  the  vector  potential  of  a  dipole  iss0 


A(i)  ~  ~.tfi  x  V' 

4* 
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(216) 


For  a  transmitter  loop  of  radius  ai,  with  unit  vector  rfi  normal  to  the  plane  of  the  loop  and  current  I\ 
flowing,  we  get 

*V>)  =  £  [a-?*]  *  v.jgijji  m 


**  Dipole  moment  is  fit,  source  coordinates  ore  primed,  field  coordinates  are  uuprimed. 
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Using  Equations  212  and  214  and  realizing  that  A  (aij)  is  approximately  constant  over  the  area  of  the 
receiver  loop,  we  find80 

Mu  =  ^—-alaln2  •  V2  x.  (nx  x  Vij-„  ■*  vj  (218) 

Three  special  cases  arise.  For  all  three,  the  distance  between  loop  centers  is  l  and  l  >  oi,oj. 

1.  Coaxial  loops  (Fig.4.5a):  We  have  •  n2  =  1  and 

Mi*  =  (219) 

2.  Coplanar  loops  (Fig.4.5b):  We  have  rij  •  n2  ~  1  and 

Mu  =  (220) 

3.  Axes  of  two  loops  intersect  at  right  angles  (Fig.4.5c):  We  have  ni  •  ri2  =  0  and 


Mu  =  0  (221) 

Even  when  the  loop  separation  is  not  many  times  the  coil  diameters,  these  general  trends  hold.  The 
mutual  inductance  of  the  coaxial  geometry  is  greater  in  magnitude  and  opposite  in  sign  to  that  of  the 
coplanar  geometry  and  the  mutual  inductance  of  the  orthogonal  axes  geometry  is  substantially  smaller 
in  magnitude  than  the  other  two.  For  this  reason,  cases  1  and  2  are  often  called  “maximum  coupled” 
while  case  3  is  called  “minimum  coupled" . 

According  to  Equation  209,  any  change  in  the  voltage  induced  in  the  receiver  coil  due  to  the  presence 
of  an  anomaly,  could  be  interpreted  as  a  change  in  the  mutual  inductance  of  the  tiansmitter/receiver 
system,  i.e., 

AV,  _  A  Mu 
V2  Afi2 

This  ratio  is  called  the  “electromagnetic  anomaly”  (measured  in  ppm)  and  it  should  be  remarked  that 
the  quantity  AMu  is  generally  complex  since  V2  in  the  presence  of  the  anomaly  is  generally  out  of  phase 
with  V2  in  the  absence  of  any  anomaly. 

Such  small  changes  in  mutual  inductance  would  stand  out  beat  for  a  system  with  a  small  mutual 
inductance  to  begin  with.  This  suggests  a  minimum  coupled  system.  This  is,  of  course,  only  true  in  air. 
If  one  wishes  to  detect  an  anomaly  (object)  buried  in  a  conducting  medium,  it  is  necessary  to  determine 
the  mutual  inductance  due  to  the  presence  of  the  medium.  This  involves  detailed  electromagnetic  analysis, 
as  we  shall  see  in  the  next  section,  and  is  very  dependent  on  geometry.  As  an  example,  if  the  medium 
is  assumed  to  be  an  infinite  flat  sheet,  the  mutual  inductance  can  be  determined  for  a  variety  of  coil 
configurations.  These  are  shown  in  Figs. 4. 0a  to  4.6e. 

Clearly  for  a  viable  system,  wc  want  a  geometry  that  maximises  the  change  in  inductance  for  the 
object  of  interest  while  minimising  that  due  to  the  host  medium.  Thus  for  example,  the  configuration 
of  Fig.4.6a  might  be  preferred  to  detect  a  small  object  buried  just  under  the  surface,  whereas  that  of 
Fig.4.6e  would  be  preferred  to  detect  the  conductive  sheet  itself.  The  choice  of  optimum  geometry  Is 
often  very  difficult  to  make  and  iu  many  cases,  it  is  unclear  which  of  several  geometries  is  optimum,  if 
indeed  any  are. 

4.2.4  Detection  Techniques 

We  have  been  considering  up  to  now  only  inductively  coupled  systems.  However,  limited  applicability  has 
been  found  for  conductive  input-inductive  output  and  inductive  input-conductive  output  systems  which 
arc  essentially  a  hybrid  of  electromagnetic  induction  and  resistivity  method*.  These  methods  are  rather 
obscure  and  we  will  say  no  more  about  them. 

oj  are  ike  unit  nonut.1  and  radius  of  Ihc  receiver  loop. 
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Techniques  of  inteiest  in  most  applications  involve  measuring  both  the  amplitude  and  phase  of  the 
receiver  signal  relative  to  the  transmitter  (for  single  frequencies).  As  we  shall  see,  the  receiver  signal 
generally  is  not  in  phase  with  the  transmitter  signal  and  thus  the  maximum  amount  of  information 
is  derived  from  such  a  measurement.  Transient  measurements  also  derive  maximum  information  since 
effectively  they  measure  phase  and  amplitude  over  a  wide  range  of  frequencies. 

In  geophysical  applications,  there  are  two  additional,  lesser-used  techniques.  Dip-angle  or  tilt-angle 
techniques  measure  the  orientation  of  the  total  (primary  plus  secondary)  field  by  finding  the  orientation  of 
the  receiver  coil  at  which  the  receive  signal  is  a  minimum.61  Also,  one  can  make  intensity  measurements 
without  phase  measurement*,  This  is  sometimes  done  by  using  a  very  large  loop  of  wire  or  a  very  long 
straight  wire  ohmicall*;  ...oupled  to  the  ground  at  both  ends  (length  ~  kilometers)  as  a  transmitter.  A 
receiver  coil  can  then  s v.rvey  a  large  area,  travelling  up  to  1/4  to  1  / 2  the  length  of  the  wire  from  the  source 
without  need  of  an  attached  transmitter.  Analysis  for  this  method  is  difficult  compared  to  a  dipole  source 
due  to  the  presence  of  conduction  currents  in  the  ground.  For  further  descriptions  of  such  techniques, 
the  reader  is  referred  to  [4]. 

Finally,  we  must  say  a  word  about  magnetotelluric  methods.  In  the  previously  described  techniques, 
the  source  fields  have  been  provided  by  artificial  means.  It  is  possible  to  use  the  geoelectromagnetic 
spectrum  as  a  source,  however,  since  we  have  seen  that  it  is  measurable  at  low  frequencies.  The  method 
involves  measuring  the  geoelectric  field  in  one  horisontal  direction,  (z),  and  the  geomagnetic  field  in  a 
horisontal  direction  perpendicular  to  the  first  (y)  at  the  air-ground  interface.  For  a  uniform  earth  and 
plane  wave  incidence  it  can  be  shown  that6* 
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where  0  is  the  angle  of  incidence  of  the  geoelectromagnetic  plane  wave  and  p,  tr,  *  are  the  electromagnetic 
parameters  of  the  earth.  Under  these  assumptions  the  surface  impedance  Z(u)  is  independent  of  0,  which 
is  highly  desirable  from  the  poiut  of  view  of  measurements.  Recalling  our  electrostatics  discussions,  we 
recognise  that  Equation  223  can  be  arranged  to  yield  a  formula  for  resistivity,  p  -  1/cr, 
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and  we  see  that  if  the  earth  is  horisontally  stratified,  we  can  deduce  information  about  it  by  measuring 
the  effective  resistivity  given  by  Equation  225.  As  for  conductivity  measurements,  analysis  can  be 
very  difficult  and  results  difficult  to  interpret  depending  on  the  complexity  of  the  underlying  strata. 
Models  have  been  developed  for  different  plane  wave  polarisations,  various  types  of  stratification  and  for 
a  current  sheet  source  at  ionospheric  altitudes.  The  main  use  of  magnetotelluric  measurements  is  for 
depth-sounding,  since  most  models  assume  that  horisontal  field  variations  are  small. 


4.3  Quasi-static  Electromagnetic  Theory 

We  start  off  as  usual  with  the  basic  Maxwell's  equations, 

(226) 

V  •  J 0=p  (227) 

*'  The  minimum  it  •  null  for  tree  iptti.  In  the  presence  of  a  conductor,  the  primer?  end  seconder?  fields  ere  out  of  phesc 
end  hence  we  get  en  elliptical!?  polerited  total  field  which  ceanot  be  represented  by  e  simple  vector.  We  gel  »  minimum 
when  the  raejor  exit  of  the  polerlsetlcn  ellipse  it  conlelned  within  the  plane  of  the  receiver  loop. 

**  It  is  estumed  thet  ea  Incident  plene  weve  hes  e  magnetic  component  in  the  y  direction  end  t revels  in  the  •  —  a  plane. 
The  earlh/slr  Interface  is  in  the  t»  -  y  plane  at  m  =  0.  See  [3]  Chapter  VI  for  details. 
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VB  =  0 

(229) 

and  we  assume  linear  isotropic  media  such  that 
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J  =  crE 

(232) 

Furthermore,  we  assume  that  charge  does  not  accumulate®3  and  hence  p  = 

0.  Maxwell’s  equations  become 

„  *  8H 

(233) 

V-E  =  0 

(234) 

V  x  H  =  aE  +  +  [jo] 

(235) 

V-H  =  0 

(236) 

The  current  density,  Jo,  in  square  brackets  in  Equation  235  is  added  when  there  is  an  impressed  current, 
that  is,  a  current  source  independent  of  the  electromagnetic  held.  Equations  233  and  235  can  be  combined 
to  yield 
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— - eu-— -  =  0 

dt  *  fit2 

(237) 

dH  d*H  n 

'fit  -</ifit2  -° 

(238) 

Clearly,  E  and  H  propagate  as  dispersive  waves.  If  we  assume  time  harmonic  solutions  for  E  and  /?,  (i.e. 
H  (f,  t)  =  !R  (f,w)cxp(iwt)|)  then  Equations  237  and  238  become 


VlE  -  itrpuM  +  e/xu»* JS?  =  0  (239) 

V*/?  —  itrfivii  +  =  9  (240) 

It  is  cumbersome  to  deal  with  two  quantities  which  satisfy  Equations  237  and  238  and  which  must 
also  satisfy  234  and  238.  It  would  be  convenient  if  E  and  U  could  be  replaced  by  a  single  potential.  There 

are  numerous  potentials  from  which  to  choose,  and  generally  we  choose  the  one  which  most  simplifies  a 

particular  problem.  For  the  present,  we  whall  choose  the  vector  potential,  >1,  given  by 

B  -  V  x  A  (241) 


which  satisfies  Equation  229.  Also  by  substituting  241  into  226  we  get 

E  =  -~-V<fi  (242) 

where  0  is  the  scalar  potential  associated  with  A.  Since  we  will  be  restricting  ourselves  to  source- free 
problems,  it  is  wise  to  choose  the  Coulomb  gauge. 

V  •  A  =:  0  (243) 

**Fur  a  finite  conductivity,  charge  dcntily  reachct  ita  •Uady-ttalc  value  in  a  time  I  ft)  t/cr. 
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(244) 


Substituting  Equ&tion  234  into  242  and  using  243  yields  <j>  =  0  and 


If  we  substitute  Equations  244  and  243  into  235,  we  find  (neglecting  J0) 

VJA  -  <rn—  -  en-Qp  =  0  (245) 

or  in  time  harmonic  form, 

V3A  -  ia-fiuA  +  e(Ui)3A  =  0  (246) 

Now  we  must  look  at  the  relative  contributions  of  the  terms  in  Equation  246.  In  free  space  po  = 
4x  x  10~7  H/m,  «o  as  8.89  x  10-1J  farads/m  and  a  «  0.  Typical  induction  frequencies  are  from  1-100 
KHi.  Thus  Equation  246  becomes 

V3A  =  0  (247) 

The  ratio  of  magnitudes  of  the  third  to  the  second  term  of  Equation  246  is  equal  to  (e/<r)w  where 
c  =  (r(o  and  e,  is  the  relative  permittivity  (dielectric  constant)  of  the  material.  A  look  at  Table  3.1 
shows  that  er  is  typically  between  1  and  10  with  the  exception  of  water  which  is  81  (rocks  typically  have 
er  ~  9,  ice  has  (,  as  4  and  er  for  soils  varies  from  4-30  depending  on  water  content).®4  Thus  even  for  an 
extreme  case  when  <r  ~  10"*  S/m  and  er  =  9,  the  ratio  of  the  third  to  second  term  is  ~  5.6  x  10-4  at 
1  KHi  frequency.  Even  at  100  KHt,  the  ratio  is  only  0.056.  For  highly  conductive  rocks  with  a  ~  104 
S/m,  the  ratio  drops  to  5  x  10"u.  Note  that  in  water  (r,  =  81, <rm„  ~  3  x  10~J  S/m)  or  some  very  wet 
soils  (efm«.  ~  30,  ~  10"1  S/m)  the  real  term  may  be  non-negligable  compared  to  the  imaginary 

term  at  the  high  end  of  the  induction  frequency  range  (~  100  KHz).  For  the  most  part,  however,  we  can 
neglect  the  third  term  in  Equation  246  and  write 

V3A  -  uryjuA  =  0  (248) 

or  in  the  time  domain, 

V*A-<^  =  0  (249) 

The  latter  equation  is  recognized  as  the  vector  diffusion  equation  and  thus  A  does  not  propagate  as 
a  steady  wave.  The  neglection  of  the  third  term  in  equation  Equation  246  is  equivalent  to  neglect¬ 
ing  displacement  current.®*  This  in  turn  is  equivalent  to  saying  that  electric  and  magnetic  fields  are 
instantaneously  propagated.  Equation  249  can  also  be  written  as 

V3A  =  -nJ  (260) 


whose  solution  in  unbounded  space  is  given  in  rectangular  coodinates  by*® 


(251) 


with  J  =  0  outside  the  volume  V'. 

44  Dielectric  constant*  vary  ilowly  in  the  XHi  to  MHs  region  tlnce  any  reionant  itructure  occur*  at  much  higher  frequencies. 

"Tbi*  approximation  i*  generally  referred  to  a*  the  “quasi -italic  approximation1*.  There  it  tome  confusion  in  the  u*e  of 
this  name.  Sometime*  the  name  refer*  to  the  situation  when  frequency  1*  to  low  that  =  0.  Alio,  let  d  he  a  typical 
source  dimension,  r  the  source  to  detector  distance  and  A  the  wavelength.  Jackson  [5]  defines  the  “near  (static)  tone”  a* 
being  d  <  r  <  A  and  the  “quasi-static  cone"  as  being  r  <  A.  However,  Walt  (3]  says  that  the  “quail-static  limit"  occur* 
in  what  Jackson  calls  the  near  aone. 

44 Equation  251  it  not  generally  valid  for  curvilinear  coordinate*  in  the  tame  manner  that  the  components  of  Vs  A  have 
no  simple  meaniug  for  other  than  cartesian  coordinates. 
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Likewise,  the  same  approximations  make  Equations  237  to  240  become 


v’l 1  =  =  0  (252) 

VSJ5  =  -urfiuE  =  0 

=  -urfiuH  =  0  (263) 

To  illustrate  the  general  nature  of  solutions  of  Equations  252  and  263,  we  choose  a  situation  in  which 
the  fields  vary  only  with  x  and  in  which  H  is  plane  polarized  in  the  y  direction.  Then 

t)  =  Hq  exp  [iwi  ±  (254) 

Jj(M)  =  <rEM{z,t)  =  (iaitu)^9  Hy{z,t)  (255) 

This  may  be  written 

Hy  =  if0exp  |-  {—■)  '  a]  exp  (wi  -  {“}  1  ®j]  (256) 

JM  =  (w7ttt)1/Jir0exp  |-  {p~- )  '  ssj  exp  j^i  (ui  -  \p7pf  '  (257) 


These  equations  represent  a  highly  damped  dispersive  wave.  For  (t/mu  1J  (l  is  the  thickness  of 
the  conductor  in  the  x  direction),  the  field  barely  penetrates  the  conductor  and  the  induced  current 
responsible  for  the  decay  of  the  field  also  is  concentrated  near  the  surface.  As  oyrw  — ►  oo  the  induced 
current  becomes  a  surface  current.  For  small  <r\iu  the  magnetic  field  penetrates  the  conductor  with  little 
attenuation  as  does  the  induced  current.  The  magnitude  of  the  latter,  however,  is  vanishingly  small. 
For  intermediate  values  of  <rpw,  a  moderately  strong  current  field  is  induced  in  the  conductor  which 
appreciably  alters  the  magnetic  field-  The  currents  are  generally  not  in  phase  with  the  magnetic  field 
and  are  concentrated  near  the  surface  of  the  conductor.  These  general  trends  are  also  found  to  hold  true 
for  three-dimensional  fields. 

The  fields,  besides  satisfying  Equations  248,  249,  262,  253,  must  also  satisfy  the  appropriate  boundary 
conditions  at  the  interface  between  different  media.  These  may  be  found  in  any  standard  text  (i.e.  [5], 
p.17)  and  are 


ft  x  (si  -  Ei)  =  0 

(258) 

ft  •  —  <TfEjj  —  0 

(259) 

ft  x  (jh,1  -  Hi )  =  K 

(260) 

ft-^ffi-pjJffj)  =0 

(261) 

n  is  the  normal  to  the  interface  between  medium  1  and  2.  if  is  a  surface  current  which  can  flow  if  either 
<T\  or  <rj  is  infinite.  Note  that  K  must  be  parallel  to  the  interface  surface  at  all  points. 

4.4  Electromagnetic  Induction  Modelling 

As  with  magnetostatics  and  electrostatics,  it  is  desirable  to  be  able  to  calculate  the  response  of  an 
electomagnetic  induction  system  to  conducting  bodies  of  arbitrary  shapes.  Unfortunately,  the  number  of 
models  which  yield  tractable  solutions  are  small,  For  instance,  for  a  dipole  source  such  solutions  can  be 
found  only  for  the  sphere,  the  thin  horizontal  sheet,  infinitely  conductive  half-plane,  and  the  horizontally 
layered  half-space.  We  will  first  choose  one  of  these,  the  sphere,  to  illustrate  the  methodology.  Afterwards, 
we  will  discuss  methods  which  can  solve  arbitrary  geometries  with  the  use  of  a  computer. 
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4.4.1  Simple  Circuit  Model 

Before  we  begin,  however,  it  is  instructive  to  investigate  an  overly  simplistic  model.  This  model,  although 
not  very  useful  for  analysing  realistic  situations,  demonstrates  most  of  the  aspects  of  electromagnetic 
induction  without  obscuring  results  in  the  mathematics  of  electromagnetic  theory.  The  geometry  is 
shown  in  Fig.  4.7.  The  method  follows  Grant  and  West  [4]. 

The  electromagnetic  induction  system  is  assumed  to  consist  of  fixed  separation  transmit  and  receive 
circuits.  The  transmitter  current  is  time-harmonic  and  both  in-phasc  and  quadrature  receiver  components 
are  measured.  The  target  is  an  isolated  conductive  circuit  with  lumped  resistance,  R,  and  inductance, 
L.  If  we  label  the  transmitter  circuit,  “0”,  the  target  circuit,  “1”,  and  the  receiver  circuit,  “2”,  then  the 
mutual  inductance  between  any  two  circuits  i,j  is  denoted  Afy  for  i,j  =  0, 1, 2.  It  should  be  noted  that 
the  difficult  aspect  from  a  field  theory  point  of  view  is  to  calculate  the  Afy. 

Assume  that  the  current  in  the  transmit  loop  is  I oe*"*.  Then  the  emf  induced  in  the  reviver  by  the 
primary  field  is 

V±p)  =  ~iuMo3I0eiut  (202) 

and  the  emf  induced  in  the  target  circuit  is 

Vi  =  -iuMoiIoe^  (263) 


The  inductance  of  the  loop  generates  a  back  emf  and  this  plus  Vi  plus  the  resistive  dro\  must  equal  zero. 
Thus 


Vt  +  Vt'  =  Vi  -  Rheiut  -  iu>Lheiut  =  0 


(264) 


Some  rearrangment  yields  the  current  in  the  target 


he 


iuL(R  —  iuL) 
R2  +  u2L2 


I0eiut 


(266) 


The  emf  induced  in  the  receiver  loop  due  to  the  secondary  magnetic  field  generated  by  the  current  in 
the  target  loop  is 

V,{,)  =  -iu>Ml3heiut  (266) 

We  define  the  “response”  as  the  ratio  of  secondary  to  primary  voltages  in  the  receiver  loop.  Using 
Equations  262,  285  and  266  and  simplifying  we  get 


<?(«)  = 


Or*  4-  ia 

1  4-  a* 


(267) 


where 


a  =  vL/il 


(268) 


and 


AfpjMn 
Afos  L 


(269) 


Studying  Equation  267  further,  we  see  that  the  coefficient  0  depends  only  cn  the  relative  sise  and 
positions  of  the  circuits  while  the  term  in  square  brackets  is  a  complex  quantity  which  is  a  function  of 
frequency  and  the  electromagnetic  properties  of  the  target. 

We  call  0  the  coupling  coefficient.  It  can  alternati  'ely  be  expressed  as 


k<n 


(270) 


where 

i  =  0,1,2  (271) 

and  |hy|  <  1.  We  see  that  each  ktj  s  the  coupling  coefficient  for  the  corresponding  circuit  ij.  Thus  the 
coupling  coefficient  0  measures  the  ratio  of  flux  coupled  into  the  receiver  via  the  target  to  that  which 
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directly  couples  to  the  receiver.  Its  value  changes  only  with  the  position  of  the  system.  The  complex 
function 

f(Q)  =  I^  =  *(Q)  +  ,T(a)  (272) 

is  called  the  “response  function”  and  a  is  the  “response  parameter”.  Since  inductance  is  related  to 
loop  diameter,  the  only  loop  parameters  which  affect  the  response  function  are  resistance  and  sise.  The 
response  function  is  shown  in  Fig.  4.8. 

Note  the  limits  of  F(a).  As  a  — »  oo,  X(a)  — » 1  and  Y(a)  -*  0.  This  is  referred  to  as  the  “inductive 
limit”  and  the  response  is  given  by 

(273) 

As  a  — ►  0,  F(a)  — »  ia.  This  is  the  “resistive  limit”  and  the  response  is 

Gr  =  ip a  (274) 

Thus  for  low  response  parameter  values,  the  response  is  small  in  amplitude  and  is  in  quadrature  with 
the  primary  field.  The  amplitude  increases,  initially  linearly,  until  the  inductive  limit  is  reached.  The 
phase  goes  from  90°  at  the  resistive  limit  to  0°  at  the  inductive  limit.  The  reasons  for  these  qualitative 
trends  are  fairly  simple  to  explain.  When  a  is  small,  little  current  is  induced  in  the  target  loop  and 
the  secondary  field  is  much  smaller  everywhere  than  the  primary.  Since  these  two  fields  do  not  strongly 
interact,  we  can  treat  each  induction  process  independently.®7  The  secondary  emf  is  shifted  180°  due  to 
two  induction  processes  whereas  the  primary  emf  is  only  shifted  90°  due  to  one.  Thus  the  response  for 
small  wL/R  is  small  and  in  quadrature  with  the  primary  emf.  As  a  increases,  the  secondary  field  induces 
an  emf  in  the  target  loop  which  becomes  comparable  to  that  induced  by  the  primary  magnetic  field.  As 
Fig.  4.9  shows,  the  phase  of  the  current  in  the  target  loop  and  thus  the  phase  of  the  secondary  field  must 
shift  so  that  the  net  induced  emf  balances  the  resistive  loss.  At  the  inductive  limit,  the  emfs  induced  in 
the  target  loop  by  the  primary  and  secondary  fields  become  equal  (see  Equation  264).  Thus  the  primary 
and  secondary  fields  must  be  in  phase  but  opposite  in  sign. 

The  inductive  limit  corresponds  to  R  — ►  0  or  infinite  conductivity.  For  a  perfectly  conductive  object, 
the  magnetic  field  cannot  penetrate  the  surface  of  the  medium.  This  happens  because  a  secondary 
magnetic  field  is  generated,  equal  and  opposite  to  the  primary  field  such  that  the  total  field  is  sero  inside 
the  object.  In  this  simple  circuit  model,  this  occurs  when  the  primary  magnetic  flux  equals  the  secondary 
flux,98  making  the  total  flux  through  the  circuit  sero  even  though  the  total  field  at  any  point  is  not 
necessarily  sero. 

Next  we  investigate  the  coupling  coefficient  0  further.  The  typical  magnitude  of  the  response  of  a 
horisontal  loop  system  passing  over  a  buried  conductor  is  shown  in  Fig.  4.10.  To  illustrate  that  our  simple 
model  predicts  such  a  response,  we  calculate  0  using  flux  diagrams,  Fig  4.11  and  4.12,  in  which  the  loop 
directions  are  chosen  such  that  hot  is  positive.  It  is  assumed  that  the  transmitter/receiver  separation 
is  fixed  so  that  hot  is  constant.  From  Fig.  4.11  we  see  that  when  the  loops  straddle  the  target,  0  (see 
Equation  270)  is  negative  and  from  Fig.  4.12,  when  the  loops  are  to  one  side  of  the  target,  0  is  positive. 
When  one  loop  is  directly  over  the  target,  0  —  0.  Thus  the  response  should  look  like  Fig.  4.10  with  each 
sero  crossing  corresponding  to  a  position  where  one  loop  is  directly  over  the  target.  The  magnitude  of  0 
is  larger  for  the  negative  portion  of  the  response  since  the  coupling  is  stronger  when  the  loops  straddle 
the  target  as  opposed  to  being  on  one  side  of  it.  Also,  the  response  is  symmetrical  about  the  target 
position,  due  to  the  fact  that  the  coupling  coefficient  is  independent  of  which  loop  carries  current. 

Finally,  if  we  express  the  response,  G(w)  (given  by  Equation  267),  in  the  Laplace  domain  and  obtain 
the  inverse  Laplace  transform  of  G(s)/s,  we  obtain 

g'{t)  =  L'x  {G(s)/s}  =  ir*  (t  >  0)  (276) 

where  r  =  L/R.  This  is  the  response  due  to  a  unit  step  current  and  it  is  a  single  damped  real  exponential 
with  time  constant  equal  to  that  of  the  target  loop. 

*rPul  mother  way,  the  back  emf  In  the  target  loop  can  be  ignored. 

MAt  the  Inductive  limit  V/  a  -iuLIitiu>.  Substitution  of  this  and  Equation  303  in  304  yield*  LI\  +  Afot  Jo  =  0. 
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The  response  characteristics  of  the  simple  circuit  model  are  similar  to  the  response  characteristics  of 
a  number  of  more  realistic  models.  The  general  shape  of  the  response  function  is  similar  to  Pig.  4.8, 
although  the  equations  governing  F(a)  are  more  complicated  than  Equation  272.  Generally  there  are 
an  infinite  number  of  poles  in  the  response  function  rather  than  one  and  this  manifests  itself  in  a  time 
response  which  is  on  infinite  sum  of  damped  exponentials  rather  than  a  single  exponential.  The  coupling 
coefficient  is,  as  mentioned,  geometry  dependent  and  must  be  examined  separately  for  each  particular 
model.  However,  for  the  borirontal  geometry  mentioned  previously  and  a  vertical  conductor,  the  profile 
(response  as  a  function  of  position)  behaves  in  a  manner  similar  to  Fig.  4.10. 


4.4.2  Response  of  a  Sphere  in  the  Field  of  e  Coil 

When  we  are  dealing  with  a  finite  body,  the  complete  electromagnetic  theory  is  necessary.  We  will  treat 
the  case  of  a  homogeneous  sphere  in  the  field  of  a  loop  as  an  example  of  such  an  analysis.  The  analysis 
has  the  advantage  of  being  a  realistic  geometry  for  detection  of  nearby  small  objects  while  reducing  in 
the  limit  to  the  commonly  used  uniform  primary  field  or  dipole  source  approximations. 

We  assume  a  spherical  polar  coordinate  system  (r,  6,  <f>)  with  origin  at  the  sphere  centre  and  will 
consider  geometries  with  axial  symmetry.  The  vector  potential  A  is  a  solution  of  Equation  249  (Section 
4.3)  and  with  the  assumption  of  axial  symmetry  must  be  of  the  form 

A  =  A{r,6)4>  (276) 

where  A  is  a  scalar  and  4>  is  the  unit  vector  in  the  tf>  direction.  We  recall  that  the  Laplacian  operator 
acting  on  a  vector  can  only  be  evaluated  component  by  component  in  a  cartesian  coordinate  system. 
Thus  we  write 

A  =  -  A  sin  <f>&  +  A  cos  <j>y  (277) 

where  i,  y  are  the  appropriate  unit  vectors.  Operating  on  each  cartesian  component  with  V3  and  com¬ 
bining  them  yields,  after  some  manipulation, 


Equation  249  thus  becomes  a  scalar  linear  2nd  order  partial  differential  equation. 

=  °  (279) 

r*  sin*  8  at 

As  usual  we  solve  this  equation  by  separation  of  variables  and  assume  time  harmonic  solutions.  Thus 


Substitution  in  Equation  279  leads  to 


A{r,8,i)  =  R(r)Q(8)eiut 

(280) 

2“* 

(281) 

(282) 

where  u  =  cos  8,  k3  =  irrpw  and  n  is  a  separation  constant. 

Equation  281  is  an  associated  Legendre  equation  with  real  solutions  /^(u),  Qa{u)'  Equation  282  is 
a  modified  Bessel  equation  with  two  solutions  (Ar)  and  v"1/*/.,,,..  j(Af)  (4^0).  If  k  =  0 

(nonconductive  region)  Equation  282  reduces  to 


45 


and  has  the  general  solution  Crn  +  Dj*-(n+1)  with  C  and  D  being  constants.  The  general  solutions  of 
Equation  279  are  then 


A(r,8,i)  =  r“1/2  [cin+^(kr)  +  D/_B_|(Ar)]  [EPi  (cos 6)  +  FQ\  (cos 6)]  (284) 


for  k  ^  0  and 


A(r,  0,t)  =  [Crn  +  [EP^  (cos  6)  +  FQ\  (cos  5)] 


(285) 


for  h  =  0.  We  can  now  turn  to  the  problem  of  the  sphere  in  the  field  of  a  loop.  The  geometry  is  shown 
in  Fig.  4.13. 

The  sphere  of  radius  a,  permeability  fi,  conductivity  <r  is  placed  in  a  nonconducting  nonmagnetic 
medium  at  the  centre  of  the  spherical  coordinate  system  ( r,8 ,  <f>).  A  circular  loop  of  wire,  radius  Rt,  is 
situated  at  (ro,0o)  and  carries  current  Ieiut.  It  is  coaxial  with  a  sensing  loop  of  radius  Rs,  situated  at 
(fs,  0s).  Both  axes  pass  through  the  sphere  centre. 

The  vector  potential  of  the  loop  is  calculated  in  the  same  manner  as  was  done  in  Section  4.2.2.  In  the 
present  case,  however,  the  current  density  is  given  by 

Jj,  =  I  sin  6' 6  (cos  8'  -  cos  do)  — —  (286) 

>•0 

Using  the  same  analysis  as  before  we  find  that  the  vector  potential,  A,  due  to  the  primary  loop  current 


is 


.«» 


.71  —  1 


$  (r<  r0) 


(287) 


Inside  the  sphere  the  vector  potential  is  from  Equation  284  70 


=  r-1/1 


£AnPt(coB0)In+i(kr) 

1.71  =  1 


$  (v<o) 


(288) 


The  secondary  or  anomalous  potential  due  to  the  presence  of  the  sphere  is  from  Equation  28571 


A<‘>  = 


•  00 

Y.EnPK  coafljr"^1) 

.71  =  1 


4>  (j*  >  o) 


(289) 


The  constants  An,  En  are  determined  from  the  boundary  conditions  (Equations  258  -  261)  of  which 
Equations  258  and  260  may  for  the  present  problem  be  rewritten, 

AWa  aW+aW 
ISAM  1  (BAM  . 


H  Sr 


-  _L  ( £d_l  + 

~  Ho  \  8r  8r  ) 


(290) 

(291) 


Substitution  and  simplification  yield  the  anomalous  or  secondary  potential  due  to  the  sphere 


=  j  Hoi  sin  Op 


*31 


;X*(ha)Pl  icoi6)  Pi  (™0o) 


where 


and  as  usual 


2n(n  +  l)rgr'*+l  ‘ 

.  .  _  ((»  +  l)/*t  +  w]  /B+t(*a)  ~  kal^^ka) 

**  °  ”  »(/v  -  + 


} 


l 


H  —  PrPO 

Ho  =  4tr  x  10"7  Il/m 


(292) 

(295) 

(294) 


Special  Cases  of  Equation  292  t 

**Tlrae  dependence  1»  implicitly  aaaumed  lo  be  and  has  been  dropped  In  (he  following  analysis. 

T0Nole  that  Q),and  f_n_  |  have  been  omiUed  since  (he  former  la  lingular  for  coal  a  1  and  (he  latter  la  lingular  for  r  =  0. 
u  Note  that  the  r*  form  of  the  solution  haa  been  dropped  since  it  la  singular  atraoj. 
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1.  Dipole  source  -  If  d  Rt  the  field  of  the  loop  should  become  dipolar.  By  expanding  (cos  #o) 
as  a  power  series  in  (1  —  cos  do),  it  is  easy  to  show  that, 


it*)  _ 

Adipolt  - 


Momr  v-'  ra2n+1 


where 


[  ; tiMn+3 

mr  —  icR^I 


Xn{ka)Ph  (cos  8) 


<f> 


(295) 

(296) 


is  the  dipole  moment  of  the  transmit  loop. 

2.  Homogeneous  Primary  Field  -  If  only  the  n  =  1  term  dominates,72  Equation  292)  becomes 


where 


A‘Lo,  =  -4^2 sin0  2xaSir°x 

![mo  (l  +  k3a3)  +  2/s]  sinh(fca)  -■  (2/r  +  Mo)  Aacosh(fca) 
[po  (1  +  k3a3)  —  n]  sinh(fco)  -f  (m  —  Mo)  fcacosh(Ita) 

ff  -IR^ 

H“  -  2~ri T 


Note  that  the  potential  is  that  due  to  a  dipole 

m-  ~2xa3R0(X  +  iy)z 

where  X  +  *Y  is  the  quantity  in  {}  in  Equation  297. 


(297) 

(298) 

(299) 


It  is  instructive  to  examine  the  resistive  and  inductive  limits  for  the  case  of  the  homogeneous  primary 
field.  By  using  series  expansions  for  cosh,  sinh,  it  can  be  shown  that  as  |feo|  -» 0, 


m*  -  r°  (-«-)  * 

which  when  M~  Mo  reduces  to 


mnW 


(300) 


(301) 


The  first  term  of  Equation  300  is  due  to  the  induced  magnetisation  of  the  sphere  and  the  second  term  is 
due  to  the  induced  eddy  currents. 

If  jieaj  -*  00  then  we  get 

mi  =  -2ira3/fo  (302) 

Both  Equations  300  and  302  can  be  obtained  by  means  which  do  not  involve  boundary  value  methods. 
This  provides  some  insight  into  the  physical  processes  involved  and  is  also  very  useful  since  often  there 
is  no  closed  form  solution  for  many  problems. 

To  obtain  Equation  300  we  assume  initially  Mr  -  1.  In  the  resistive  limit,  only  the  primary  field 
contributes  to  the  induction  process.  We  divide  the  sphere  into  annular  rings  of  radius  rsinfl  and  cross- 
sectional  area  rd&dr.  The  emf  induced  in  each  ring  is 


V  —  -iwpoirr2  sin2  0Ho 

We  equate  this  to  the  resistive  drop  in  the  ring  and  get 

d2/  =  -^wr/iowr2  sin  9drd$Ho 

ft 


(303) 


(304) 


T2Ihl»  implies  lhal  a  <  (itj,  +  <P)X^  or  a  <  (flj  +  (l2)'73. 
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(305) 


The  contribution  to  the  dipole  moment  is 

=  x  (r  sin  6) 3  d3I 

Integrating  Equation  305  we  get  301,  which  is,  of  course,  300  with  Hr  —  1.  To  get  300  for  ft,  ^  1  we 
replace  HoHq  with  the  actual  flux  density  inside  the  sphere  and  then  add  a  term^to  account  for  induced 
magnetisation.  Recall  (Section  2.8)  that  if  M  is  the  magnetisation  of  a  sphere,  Hq  the  external  field  and 
Hi  the  internal  field,  then 

Hi  -  Ho  -  1/3 M  (300) 

But 

and  so 


M  =  - — —Hi 
Mo 


fiHo 


The  dipole  moment  due  to  magnetisation  is 

m  =  ^xa3M 
3 

If  we  add  Equation  309  to  301  and  replace  Mo#o  by  of  Equation  308,  we  get  the  general  resistive  limit 
Equation  300  for  any  fi. 

The  inductive  limit  is  obtained  by  finding  the  effective  magnetisation  that  makes  the  flux  density 
vanish  everywhere  inside  the  sphere.  We  have 


(307) 

(308) 

(309) 


0  =  Bi  =  Ho  [ili  +  MlJ  =  Mo  (^Hq  +  -  Af/,^ 


or 


Ml  =  -\h0 


(310) 

(311) 


Substituting  A?$  for  M  in  Equation  309  yields  302. 

Now  let  us  examine  the  secondary  potential  of  the  sphere  in  the  field  of  a  loop  (Equation  292)  in  more 
detail.  First  we  note  that  it  is  the  voltage  induced  in  the  receive  loop  that  we  measure  rather  than  the 
vector  potential.  This  is  given  by 


V(')  =  =  2viuRs 


(312) 


where  the  circular  integral  is  performed  around  the  receiver  loop.  In  terms  of  distances  rather  than  angles 
this  becomes, 

=  2x»Mofu>- ------- -rx 

(4  + *10* 


^  am+i  Pn1(dr/[4  +  4]*)^(^/[4  +  *l]*)/v 
h  2n(n+1)  (4  +  *f )?  (4  +  4)**1  * +  *  * 

where  dx  is  d  for  the  transmit  loop  (see  Fig.  4.13),  ds  is  d  for  the  receive  loop, 

XJka)  =  Xn(ka)  +  iY%(ka) 


) 


(313) 


and  XK,  Yn  are  real  functions. 

As  in  the  simple  loop  target  example,  the  expression  ^(w)  +  »F„(w)  contains  all  the  frequency 
dependence  while  the  remainder  of  Equation  313  contains  the  geometry  factors.  The  present  response  is 
for  more  complicated  than  that  for  the  simple  case  because  we  now  have  a  complicated  sum  of  products 
of  geometry  and  frequency  dependent  terms.  The  response  parameter  is  now  <rnua?  rather  than  wL/R. 
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Examining  Equation  313  we  see  that  the  higher  order  terms  decrease  more  rapidly  with  d  than  the 
lower  order  terms73.  Thus  if  the  sphere  is  far  from  the  loop  compared  to  its  sise  (or  much  smaller  than 
the  loop  diameter),  the  secondary  voltage  is  that  due  to  a  sphere  in  a  uniform  field.  In  this  case  we  recall 
that  the  secondary  field  is  that  of  a  simple  magnetic  dipole.  As  the  loop  approaches  the  sphere,  higher 
order  terms  become  significant  and  the  secondary  field  of  the  sphere  receives  contributions  from  higher 
order  multipoles. 

The  four  lowest  order  multipole  response  functions  are  shown  in  Fig.  4,14.  Note  that  the  higher 
multipoles  reach  the  inductive  limit  at  higher  values  of  o-jwa1.  This  implies  that  the  higher  multipoles 
are  more  difficult  to  excite. 

The  time  domain  response  may  obtained  in  the  same  manner  as  for  the  loop  target  (Section  4.4.1). 
For  a  unit  step  primary  current  we  find  for  t  >  0,7i 


V(‘)(t) 


2xRtRsI 

a 


£ 


(o,  dr,  ds ,  Rt,  Rs)  ^  Anm  (p,)  exp 

m=l 


where  6nm  are  positive,  real  solutions  of 


(314) 


n  (f3*  (^nm)  +  ^nmj'n-1  (d*m)  —  0  (315) 

jn  are  the  spherical  Bessel  functions, 


(2n  +  l)a’"-^  (dr/  [4  +  Jg]  »)  ^  (ds/  [d»  +  *»]  ») 

”  n(n  +  l)(4  +  ^.)afi(d|  +  f27)^ 


(316) 


and 

A  —  _ _  ^nn»  _  .  . 

""  “  n/i,  (np*  +  1)  +  *»m  -  n(n  +  1)  (31') 

A  number  of  features  are  worth  comparing  with  the  simple  loop  target  problem.  First,  Equation  314 
is  a  sum  of  damped  real  exponentials  as  opposed  to  a  single  damped  exponential  as  in  Equation  275.  For 

the  simple  loop  target,  the  geometry  or  coupling  was  contained  in  the  coefficient  /3,  whereas  for  this  case 

the  geometry  determines  the  coefficients  RT,  Rs,  Wn.  For  the  simple  loop  target,  information  about  the 
object  sue  and  shape  was  contained  in  the  coefficient  of  time  in  the  exponential.  This  is  also  true  for  the 
sphere  in  the  field  of  a  loop  except  that  there  are  now  an  infinite  number  of  coefficients,  anm, 


- 


pKTtt5 


(318) 


In  spite  of  the  differences,  however,  we  shall  see  that  the  sphere  and  simple  loop  target  have  timo  responses 
which  appear  quite  similar,  just  as  the  frequency  responses  are  similar. 

It  can  be  shown  that  for  constant,  n,  the  increase  as  m  increases  (Table  4.1).  Thus  higher  order 
terms  for  each  multipole  decay  faster  than  the  lower  order.  For  »  =  1,  p,  =  1,  we  have  the  special  case 
of  a  noopermcable  sphere  in  a  homogeneous  field,  In  this  case 


list 


_ 4&s 

(4  +  *W4(4  +  *i)* 


iH-(s£)‘]  i‘>o)  <si9> 

and  the  initial  amplitudes  of  all  exponentials  are  the  same  (Fig.  4.15).  This  allows  the  time  constants  to 
be  extracted,  in  principle. 

For  n  —  1,  Hr  /  1,  the  initial  amplitudes  of  the  exponentials  are  initially  small  for  m  =  1  and 
increase  slowly  with  m,  asymptotically  approaching  the  value  of  the  nonpermcable  amplitudes  for  large 
m  (Fig.4.16).  This  makes  extraction  of  the  time  constants  by  simple  means  impossible.  We  will  return 
to  the  problem  of  extracting  time  constants  later. 

nTW»  U  culctl  It-  t«  If  we  Ut  d  =  i*  a  dv  i.v.  Uu>  cooU/.-f 
M  A  dtuUcd  dtriv»*ioe  t*  found  ln(6dj. 
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4.4.3  Solution  Using  Impedance  Boundary  Conditions 

One  drawback  to  boundary  value  problems,  such  as  the  previous  example,  is  that  there  are  only  a  limited 
number  of  solvable  geometries.  Numerical  methods  would  have  a  much  broader  range  of  applicability 
since  many  different  geometries  could  be  used  as  an  input.  If  the  object  to  be  detected  has  sufficiently 
high  conductivity  so  that  impedance  boundary  conditions  can  be  applied,  then  the  response  will  depend 
only  on  surface  properties  of  the  object.  It  is  well  known  that  in  the  numerical  solution  of  scattering 
problems,  the  surface  integral  equations  are  much  simpler  to  handle  than  the  volume  integral  equations. 
It  cun  be  shown  ([5]  p.338)  that  at  the  surface  of  a  good  but  not  perfect  conductor,  the  magnetic  held, 
H,  and  electric  field,  E,  are  related  by  the  impedance  boundary  conditions 

=  (320) 

where  n  is  the  unit  normal  to  the  surface  and  Zs  is  the  surface  impedance  of  the  conductor  with  electro¬ 
magnetic  parameters  «r,  p,  e. 

Zs  =  Ze=  (-^-V  (321) 

\<X  +  W€/ 

where  Zc  is  the  intrinsic  impedance.  This  is  strictly  valid  only  if  spatial  variations  of  the  fields  orthogonal 
to  the  surface  normal  nre  small  inside  the  conductor.  For  closed  bodies  this  condition  is  true  if 

tei)1'’4'’*1  <m> 

where  al!  quantities  arc  calculated  inside  the  conductor,  p  is  the  smallest  radius  of  curvature  at  the  point 
on  the  surface  and 

(323) 

is  the  propagation  constant  in  free  space.  Equation  320  may  be  written  os 

M  =  (324) 


where  Af ,  j  arc  the  magnetic  and  electric  surface  currents  which  arc  generated. 


A?  «  -  A  x  £  (325) 

/=  A  x  H  (326) 

The  general  scattering  geometry  is  shown  in  Fk.4.17.  The  incident  (primary)  electromagnetic  fields 
E*,  H'  induce  currents  an  the  surface,  S,  of  an  imperfectly  conducting  body  such  that  the  impedance 
boundary  conditions  on  5  are  satisfied.  The  scattered  fields  at  r  due  to  the  surface  currents  are  given  by 
([3]  pp.21-24) 

i-Vx  Vx.f-Vx  P  (327) 

iw<0 

//*  I  r~  v  x  V  x  f'  +  V  x  A  (328) 

wupo 

where  A  and  F  are  the  magnetic  and  electric  vector  potentials  or  Schdkunoff  potentials  given  by 


A-  Js  J  (r')  G  (r,  /)  d* 


(329) 


F=  J  Af  (r')  <?  (r,  r')  d* 
where  G  is  the  free  space  Green’s  function, 


(330) 


(331) 


SO 


Substitution  of  Equations  328,326,329,  330  into  327,328  and  recognizing  that  E,  H  in  Equations  325,326 
are  the  total  fields  given  by  E1  +  E‘,  H{  +  H‘  respectively,  we  get  two  integral  equations  [67]  : 


J{r)-2h(f)  x  H*  (f)  +  2n 


x  V'Gds' 


and 


-2 n(f)  x  jy  iu>e0M  (?)  G  -  - -J-V's  •  M  (?')  V'Gdn1 

#  (,-J  =  *2*  ^ /(,-■)  Gds'  -  ±Jt  Vi  •  f{?)  VGd.' 
-  J  M  (p)  x  V'Gdt' 


(332) 


(333) 


where  V'  is  the  tangential  gradient  with  respect  to  r'.  The  coupling  between  the  Magnetic  Field  Inte¬ 
gral  Equation  (332)  and  the  Electric  Field  Integral  Equation  (333)  is  resolved  by  using  the  generalized 
impedance  boundary  condition, 


M  (?)  =  Zs  (?)  [/  (?j  x  h  (p)]  (334) 

where  now  Zg  is  a  dyadic  matrix  (2  x  2)  and  J  is  the  surface  current  (  a  2-vector  with  components 
usually  chosen  along  the  directions  of  principle  curvature).  If  we  write, 

-w-(4o 

then  there  are  a  number  of  approximations  for  the  surface  impedance.  Two  commonly  used  ones  are: 

1.  Zg  constant  everywhere  on  5  (Leontowich  approximation). 

Zx  =  -Zt  =  Ze  (330) 

This  approximation  is  useful  when  Equation  322  is  fully  satisfied. 

2.  Zg  is  curvature-dependent. 

Zi  =  [l  -  Q  (P)]  Ze 

Z3  =  [l  4-  Q  (p)]  Zc 

Q(?)  =^(1-x)6{Km-K.)  (337) 

where  fftt,  K*  are  the  principle  curvatures  at  ?  and  u.  D  are  two  orthonormal  tangential  vectors  at 
r'  for  which  u  x  v  -  n  is  true.  This  is  useful  when  the  radius  of  curvature  on  some  portions  of  the 
body  approaches  the  electromagnetic  wavelength. 

Equations  332  -  334  sing  the  appropriate  r.pproximation  for  Zg  may  be  solved  by  one  of  a  number  of 
numerical  methods.  One  popular  one  is  thj  Moment  Method  in  which  the  integral  equation  is  solved 
numerically  by  expanding  the  incident  fields  on  the  surface  in  a  series  of  suitable  basis  functions,  applying 
a  set  of  testing  functions  and  inverting  the  matrix  of  coefficients  of  the  resulting  set  of  simultaneous  linear 
equations.  There  are  a  host  of  other  numerical  methods  that  may  be  used,  but  a  full  discussion  is  beyond 
the  scope  of  this  course.  Further  details  may  be  found  in  [69]. 


*  (-')  - 
0  / 


(335) 
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4.5  Inverse  Electromagnetic  Induction  Problems 

The  object  of  inverse  electromagnetic  induction  problems  is  similar  to  that  of  Section  2.10  for  the  magne¬ 
tostatic  case,  namely  find  the  position  (including  depth),  shape  and  size  parameters  of  a  hidden  compact 
object  of  finite  sixe.75  As  for  the  magnetostatic  problem,  we  are  faced  with  the  problem  of  uniqueness. 
Although  generally  speaking,  many  sources  may  yield  the  same  set  of  field  values,  often  in  practice  the 
geometry  is  sufficiently  well-known  a  priori  or  the  sufficient  information  is  available  about  the  source  (i.e, 
there  are  only  a  finite  number  of  sixes  and  shapes)  that  the  inverse  problem  is  uniquely  solvable.  In  any 
event,  the  success  or  failure  of  the  inverse  method,  be  it  model  fitting  or  pattern  recognition,  is  essentially 
a  test  of  uniquess. 

There  is,  however,  a  notable  difference  between  the  magnetostatic  and  electromagnetic  induction 
inverse  problems.  In  the  magnetostatic  case,  we  saw  that  the  permeability  had  little  effect  on  the 
measured  responses.  For  e.  m.  induction,  the  permeability  can  dramatically  affect  the  response  as 
can  the  conductivity. 

We  will  restrict  our  discussion  to  the  time  domain,  although  most  of  the  arguments  related  to  local¬ 
isation  apply  to  the  frequency  domain  as  well.  Furthermore,  we  shall  consider  the  detection  of  a  buried 
sphere  using  a  set  of  horisontal  coils  as  in  Fig.4.13.  The  time  domain  response  is  given  by  Equation  314, 
which  is  a  sum  of  damped  real  exponentials.  Compact  bodies  generally  give  such  a  response  and  so  the 
following  arguments  will  also  be  generally  applicable  to  compact  objects  other  than  spheres. 

Localisation  of  the  sphere  in  a  horizontal  plane  is  not  as  difficult  as  for  localisation  of  a  magnetostatic 
dipole  in  a  plane.  The  RMS  secondary  voltage  induced  in  the  receive  coil  as  it  moves  over  the  object  in  a 
horisontal  plane  can  be  thought  of  as  a  spatial  signature  of  the  object.  This  signature  is  generally  much 
narrower  than  the  magnetostatic  signature  for  the  same  geometry.  This  is  because  e.m.  induction  is  an 
active  detection  method  whose  signal  decreases  very  roughly  as  ~  r~8,  whereas  magnetostatic  detection 
is  a  passive  method  with  a  signal  fall-off  ~  r~*.  Furthermore,  in  e.m.  induction,  the  primary  field  is 
vertical  and  hence  the  magnetic  moment  induced  in  the  sphere  is  vertical.  The  effect  of  this  is  that  the 
secondary  voltage  is  a  maximum  when  the  sphere  center  is  situated  on  the  coil  axis;  i.e.  when  the  sphere 
is  directly  under  the  coils.  Foi  non-sphericol  objects,  the  maximum  secondary  voltage  still  tends  to  occur 
when  the  object  is  under  the  coils,  although  the  geometric  center  will  not  generally  lie  exactly  on  the 
coil  axis.  Typically,  it  is  possible  to  localise  an  object  in  she  horisontal  plane  to  within  one  radius  of  the 
receive  coil. 

Depth  determination  is  slightly  more  complicated.  If  we  consider  the  geometry  of  Fig.4.13  and  assume 
that  the  object  is  small  enough  w.r.t.  the  coils  or  far  enough  from  them  such  that  only  the  n  =  1  term  of 
Equation  314  contributes,79  then  we  fmd  that  the  ratio  of  voltages  induced  in  2  copianar,  coaxial  receive 
coils  of  radius  Hi  and  Ri  at  t  =  is 


r  _  (<<)  _  #(<*»  + Bj)a/> 

Vi(U)  "  R§(d»  +  K»)S/J 

which  assumes  the  same  number  of  tviac  for  each  coil.  We  see  that  in  this  simple  case,  the  voltage  ratio 
is  a  function  only  of  the  depth,  d,  of  the  object  and  not  of  the  object  parameter*  (p,  tr,  a).  Inversion  of 
Equation  338  by  some  method  allows  us  a  means  of  determining  depth.  It  turns  out  that  the  assumption 
of  a  spherical  object  and  uniform  field  are  not  that  restrictive.  Good  results  are  obtain  d  from  the  method 
even  when  the  suhere  depth  is  only  0,3  times  the  radius  ol  the  larger  coil.77  Situations,  in  which  the 
object  has  dimensions  similar  to  the  coil  diameter  and  its  distance  from  the  coil  plane  is  less  than  a  coil 
diameter,  are  typical  of  ordnance  and  archaeological  artifact  detection.  The  response  for  such  geometries 
would  be  expected  to  have  a  relatively  high  contribution  of  higher  order  multipole  fields.  However,  the 
good  depth  estimates  that  arc  obtained  show  that  the  contribution  is  small.  In  addition,  good  results  are 
obtained  for  metallic  spheroids  with  sises  similar  to  the  spheres,  independent  of  object  orientation.  Nor 
must  the  object's  geometrical  center  lie  on  the  coil  axis.  Good  results  are  usually  obtained  for  spheres  if 
the  object  ceutcr  is  displaced  horizontally  from  the  coil  axis  by  up  to  a  fraction  (typically  ~  1/2)  of  a  coil 

’'Since  v«  ere  deeliog  with  compact  objects,  we  exclude  from  thii  ditcuailon  geophyiicei  problem)  related  to  detection  of 
tend-infudte  media  «uch  u  layered  helf->pecc)  or  vertical  dyke).  Three  ere  mo«t  often  etulyeed  by  e  modelling  epproeeb 
T*i.e,  the  primer/  field  1$  uniform 
”Tlu)  tnumn  ft,  =  3fiJt  0.3Jt,  <  a  <  O.Sfij. 
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radius.  If  the  center  is  off-axis,  the  depth,  d,  determined  from  Equation  338  approximates  (d*  +  z2)1' ! 
where  do  is  the  true  vertical  depth  and  x  is  the  horizontal  displacement  of  the  object  center  from  the 
axis.  A  more  detailed  discussion  of  this  method  is  found  in  [74].  Experiments  must  still  be  done  to  verify 
if  the  method  is  applicable  to  objects  other  than  spheroids,  for  similar  geometries. 

Assuming  we  now  know  the  location  of  th'  object,  we  are  left  with  the  problem  of  identification, 
namely  determining  size,  shape  and  material  properties  (ft,  a).  As  for  magnetostatic  detection,  there  are 
two  general  methods:  model  fitting  and  pattern  classification. 

Model  fitting  includes  least  squares  fitting  and  Prony’s  method  [78].  Of  course,  model  fitting  relies 
strongly  on  there  being  an  applicable  model  and  this  is  what  often  limits  its  use.  For  a  sphere,  the 
secondary  voltage  can  be  fitted  to  the  function 

V‘(i)  =  f;  Wn(a)  £;  A„m  exp  |  —  t}  (339) 

by  a  nonlinear  method  with  o,  a  being  the  free  parameters.  Of  course  Anm,  Knm  are  functions  of  ft  but 
for  nonpermeable  objects  they  are  one  set  of  constants  and  for  fir  ~  100  or  greater,  they  are  another 
set  of  constants.7*  For  nonpermeable  objects  only  a  few  terms  are  usually  needed  in  the  summation  but 
for  permeable  objects,  hundreds  of  terms  are  needed  to  approximate  Vf*)  accurately.  There  are  also  the 
standard  problems  of  nonconvergence  and  need  for  initial  parameter  estimates  which  all  nonlinear  fitting 
methods  have  in  common.  Most  importantly,  of  course,  the  model  applies  only  to  a  sphere. 

However,  it  is  found  that  compact  conductive  objects  generally  produce  a  secondary  voltage  of  the 
form 

00 

V"(f)  =  X>ne~«»‘  (340) 

n=l 

and  thus  a  nonlinear  fit  of  the  secondary  voltage  using  a  truncated  series  like  that  of  Equation  340 
(n  =  1, 2 ,...N)  would  yield  coefficients  A*,  a*  which  would  presumably  be  functions  of  the  size,  shape, 
<r,  fi  (and  unfortunately  the  orientation,  as  well)  of  the  object.  The  problem  in  such  an  approach  turns 
out  to  be  the  fact  that  real  damped  exponentials  are  highly  correlated  functions.  Thus  the  estimates  Am, 
o'*  are  strongly  dependant  on  the  number  of  terms  ( N )  used  in  the  model. 

Another  technique  that  has  been  used  with  tome  degree  of  success  in  underwater  acoustics  and 
electromagnetics  is  Prony’s  method.  In  those  applications,  the  waveforms  involved  are  a  sum  of  damped 
complex  exponentials.  Prony's  method  is  basically  a  nonlinear  fitting  procedure  where  the  nonlinearity  of 
the  system  it  concentrated  in  a  single  polynomial  and  thus  moat  of  the  above  comments  about  nonlinear 
fitting  apply.  One  advantage,  however,  is  that  an  initial  guess  of  the  parameters  is  not  required.  As 
with  regular  nonlinear  fitting,  the  number  of  exponentials  in  the  model  has  a  very  strong  effect  on  the 
parameter  extimatea  obtained.  A  further  disadvantage  is  that  the  estimation  procedure  is  sensitive  to 
the  sampling  interval  chosen  and  very  sensitive  to  the  S/N  ratio. 

One  might  be  tempted  to  try  a  simpler  approach  which  is  often  used  in  nuclear  physics.  The  tail  of 
the  waveform  is  fitted  to  a  single  exponential,  with  the  assumption  being  made  that  the  faster  decaying 
terms  are  negligable  at  sufficiently  great  times.  The  smallest  time  constant,  cf*,  is  obtained  and  the 
estimated  term  Aie-**1*  is  subtracted  from  the  waveform  to  obtain  a  reduced  waveform.  The  tail  of  the 
reduced  waveform  is  fitted  to  a  single  exponential  and  the  process  is  repeated.  This  method  might,  in 
fact,  work  reasonably  well  for  nonpermeable  spheres  in  which  only  a  single  exponential  is  significant  for 
large  times  (Fig.  4.15).  For  permeable  spheres,  several  exponentials  contribute  at  any  time  (Fig.  4.16) 
and  the  method  will  fail.79  In  general,  objects  of  arbitrary  shape  do  not  have  waveforms  like  Fig.  4.15. 

Even  if  filling  the  waveforms  to  a  sum  of  damped  real  exponentials  did  yield  accurate  coefficients, 
relating  those  coefficients  to  the  properties  of  the  target  objert  would  be  a  formidable  task  without  a 
priori  knowledge  of  the  object's  shape  classification  (i.e.  spheroid,  cylinder,  hoi’  'W,  solid,  etc.).  Indeed, 
even  this  knowledge  would  allow  an  analytical  solution  for  the  object's  parameters  in  only  the  few  cases 
such  as  the  sphere  or  spheroid  for  which  analytical  models  exist. 

Taln  all  the  previous  modelling,  il  ha*  beta  assumed  that  ftr  1*  constant,  although  ft  U  obviously  governed  by  some 
hyitereus  curve.  However,  ia  the  earth'*  field  fit  is  typically  between  100  and  1000  for  ferrous  object*  and  hence  the  effect 
is  at  though  fit  were  constant  (see  Table  4.1). 

n  Because  of  the  correlated  nature  of  real  exponentials,  one  can  always  fit  the  tail  to  a  tingle  exponential.  However,  the 
estimated  coefficient  could  not  be  simply  related  to  p,  tr,  «. 
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If  the  numbei  of  objects  under  consideration  is  finite,  we  may  consider  pattern  recognition.  A  method 
which  has  been  successfully  used  to  classify  metallic  spheroids  [79]  will  be  described  as  an  example.  As 
previously  stated,  once  we  have  localized  the  object,  the  coil  axes  are  approximately  over  its  geometric 
center.  Referring  to  Fig.  4.13,  where  now  the  spherical  objec;  is  replaced  by  an  arbitrary  object  with  axial 
symmetry,  it  is  obvious  that  the  response  will  be  independent  of  the  azimuthal  angle  of  the  symmetry 
axis  of  the  object. 

Feature  vectors  are  obtained  from  the  electromagnetic  induction  response  by  one  of  a  variety  of 
schemes.  This  usually  consists  of  first  normalizing  the  time  response  relative  to  a  fixed  time  after  the 
falling  edge  of  the  transmitter  pulse;  segmenting  the  normalised  response  and  fitting  each  of  the  finite 
number  of  segments  to  a  simple  function  such  as  a  straight  line,  single  exponential  or  a  constant.  The 
segments  do  not  have  to  be  equal  time  intervals  but  may  be  adjusted  for  constant  energy  (area).  There 
are  advantages  and  disadvantages  for  each  feature  type  which  are  beyond  the  scope  of  this  course  (see 
[79]  for  a  discussion).  A  simple,  but  quite  successful  choice  [79]  has  been  the  mean  response  lor  each 
segment  with  equal  time  intervals.  A  feature  vector  for  a  response  then  consists  of  the  sequence  of  fit 
parameters  for  all  segments. 

We  can  use  a  continuous  parameter  classifier  very  similar  to  that  of  Section  2.10,  except  that  now  the 
manifold  is  a  function  of  only  one  continuous  parameter  8,  the  polar  angle  of  the  object’s  symmetry  axis. 
Hence,  the  cIms  prototype  for  an  object  is  a  one  dimensional  manifold  (i.e.,  a  curve)  in  an  N  dimensional 
feature  space.*50  An  example  is  shown  in  Fig.4.18  for  an  N  =  2  dimensional  feature  space. 

The  approach  is  to  approximate  the  curve  by  a  sequence  of  line  segments  which  connect  points  on  the 
curve.  This  reduces  the  problem  to  finding  the  minimum  distance  from  the  test  vector  to  a  sequence  of 
line  segments. 

The  prototype  for  class  t  and  given  8  is  a  point  defined  by  the  head  of  a  vector  denoted  |n»i  (6)),  The 
prototype  is  approximated  by  a  finite  number  of  line  segments  joining  points  |m^)  where  the  subscript 
j  indicates  that  the  prototype  feature  vector  is  evaluated  at  the  discrete  value  of  the  parameter  8  =  8j. 
We  define  a  unit  vector  \cuj)  along  the  line  segment  joining  [m^)  and  |nqj-n), 

lau)  =  l«u)/(«ul«u)*  (341) 

where 

l*u)sl«*j+i)-l«u)-  (342) 

Let  |a)  be  a  test  vector  and  define 

l#u)  =  I*)  ~  K.i)  •  (343) 

If  we  define  |y{^  to  be  the  projection  of  |y^)  onto  the  line  segment,  and  |djj)  to  be  the  vector  which 
is  normal  to  the  line  segment  and  which  passes  through  |a),  then 


l<*u)  -  hfij)  -  |v$) 

(344) 

y\'J)  =  (yul  «u)  l“u)  =  flu  l«u) 

(346) 

where 


flu  = 


Kil«u)‘ 


(346) 


If  |z)  is  a  sample  from  the  class  i  corresponding  to  the  region  of  the  prototype  curve  for  which  8j  <  8  < 
0j+\,  then  an  estimate,  8,  of  the  continuous  parameter  associated  with  |a>)  may  be  obtained  from 


0  =  +  flu  (fy+i  ~  0j) 

a  13  la  lypical. 


(347) 
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The  minimum  distance,  dij,  from  the  test  vector  to  the  region  of  the  curve  between  |  mij)  and 
is  approximated  by 


did  ~  (dij\  dij}*  (348) 

Since  the  line  segment  only  approximates  the  manifold  between  the  two  endpoints  of  the  line  segment,  the 
previous  distance  equations  are  really  only  true  provided  lies  between  these  points.  Examination 

of  Equations  345  and  346  reveals  that  this  occurs  for  0  <  qij  <  1.  If  this  condition  is  not  satisfied  then 
dtj  is  replaced  by  the  distance  from  the  test  vector  to  the  nearest  endpoint  of  the  segment.  For  speed 
of  implementation,  qij  is  calculated  first  and  based  on  its  value,  the  appropriate  calculation  of  dij  is 
carried  out. 

Finally,  the  minimum  distance,  <U,  from  the  test  vector  to  the  curve  is  approximated  by 

di  =  nun^j}.  (349) 

The  test  vector  is  assigned  to  the  class  i  for  which  >li  is  a  minimum. 

Probability  of  correct  classification  for  a  restricted  case  (responses  of  a  set  of  four  steel  spheroids 
obtained  in  a  nonmetallic  laboratory)  are  between  90%  and  97%.  One  problem  is  that  the  responses  vary 
slowly  with  depth  and  hence  to  improve  classification  success,  the  design  set  must  include  various  depths 
or  else  a  correction  for  depth  must  be  applied  to  each  response.  The  former  would  then  require  the  use  of 
a  two  dimensional  continuous  parameter  classifier  such  as  was  used  in  Section  2.10.  The  student  should 
consult  [79]  for  more  details. 

5  Conclusion 

The  detection  of  compact  objects  by  low  frequency  electromagnetics  has  been  discussed.  In  particu¬ 
lar,  magnetostatic,  electrostatic  and  electromagnetic  induction  techniques  have  been  described  with  an 
emphasis  on  locating  and  identifying  objects  rather  than  merely  detecting  them. 

Although  the  techniques  for  compact  object  detection  are  broadly  related  to  their  counterpart  geo¬ 
exploration  methods,  it  must  be  remembered  that  the  two  problems  are  quite  different  and  methods 
cannot  be  directly  carried  from  one  problem  to  the  other.  Sources  and  geometries  are  usually  radically 
different  for  both  problems  and  consequently  detectors  and  signal  processing  techniques  are  also  different. 
Geoexploration  techniques  are  covered  in  a  number  of  standard  textbooks. 

Finally,  it  is  sometimes  stated  that  low  frequency  electromagnetic  detection  of  compact  objects  is 
a  “mature  science".  !  hope  that  by  the  time  this  part  of  the  text  is  reached,  it  will  be  clear  to  the 
reader  that  there  is  st  U  much  fundamental  research  to  be  done  before  objects  can  be  reliably  located  and 
identified. 


55 


References 


LOW  FREQUENCY  ELECTROMAGNETIC  REMOTE  SENSING 


GENERAL 


Geo-electromagnetic  Theory 

[1]  -,  IEEE  Trans  GeoSci.  and  Remote  Sens.  Special  Issue  on  Electromagnetic  Methods  in  Applied 
Geophysics  GE-22,  pp.3-90,  1984. 

[2]  J.R.Wait,  Geo-Electromagnetism.  New  York:Aeademic  Press,  1982. 

[3]  Mining  Geophysics  Volume  2  -  Theory.  Tulsa  Oklahoma:  The  Society  of  Exploration  Geophysi¬ 
cists,  1967. 

[4]  F.Grant  and  G.West,  Interpretation  Theory  in  Applied  Geophysics.  New  York:McGraw-Hill,  1905. 

Electromagnetic  Theory 

[5]  J. Jackson,  Classical  Electrodynamics.  New  York:  Wiley,  1975. 

[0]  W.Panofsky  and  M.Phillips,  Classical  Electricity  and  Magnetism  Reading,  Mass.:  Addison-Wcsley, 
1902. 

[7]  W.Smyihe,  Static  and  Dynamic  Electricity  Third  Edition,  New  York:  McGraw-Hill,  1968. 

[8]  J.Stratton,  Electromagnetic  Theory,  New  York:McGraw-Hill,  1041. 

Mathematical  Techniques 

[9]  G.Arften,  Mathematical  Methods  for  Physicists ,  New  York:  Academic  Press,  1970. 

[10]  M.Abramowiti  and  I.Stegun,  Eds.,  Handbook  of  Mathematical  Functions,  New  York:  Dover,  1906. 

[11]  O.C.Zienkiewics,  The  Finite  Element  Method  in  Engineering  Science,  London:  McGraw-Hill,  1971. 

[12]  M, Chari  and  P.Silvester  ,  Eds.,  Finite  Elements  in  Electrical  and  Magnetic  Field  Problems,  New 
York:  Wiley,  1980. 

[13]  J.Tou  and  R.Gonsolei,  Pattern  Recognition  Principles.  Reading,  Mass.:  Addison  Wesley,  1974. 

[14]  J.W.Samraon  Jr.,  “Interactive  pattern  analysis  and  classification",  IEEE  Ihins.  Computers,  vol. 
C-19,  p.594,  1070. 

[15]  J.ll.Ricc,  “Experiments  on  Gram-Sckmidt  orthonormolisation",  Math.  Comp.,  vol.29,  p.326,  I960. 

Remote  Sousing 

[16]  R.N. Colwell,  Editor-in-Chief,  Manual  of  Remote  Sensing,  Second  Edition ,  Volumes  1  and  2.  Falla 
Church  Va.:  American  Society  of  Photograminetry,  1963. 


50 


MAGNETOSTATIC  METHODS 


General  Magnetism 

[17]  R.Bosworth,  Ferromagnetism.  New  York: Van  Nostrand,  1951. 

[18]  G. Garland,  Introduction  to  Geophysics.  Toronto:  W.B. Saunders,  1971  (Chapters  15  -  21). 

Theory 

[19]  J.Simkin,  C.Trowbridge,  “Three-dimensional  nonlinear  electromagnetic  field  computations  using 
scalar  potentials",  IEE  Proc  127B,  pp.368-374,  1980. 

Applications 

[20]  J.Scarsello  and  G. Usher,  “A  low  power  magnetometer  for  vehicle  detection "  ,IEEE  Trans  Mag, 
MAG- IS,  pp.  1101-3,  1977. 

[21]  G.Ioannidis,  “Identification  of  a  ship  or  submarine  from  its  magnetic  signature”, IEEE  Trans  AES, 
AES-13,  pp.327-329,  1977. 

[22]  J.McFee  and  Y.Das,  “The  detection  of  buried  explosive  objects”,  C.J. Remote  Sensing  6,  pp.104-121, 
1980. 

[23]  J.E.McFee  and  Y.Das,  “Development  of  an  intelligent  system  to  locate  and  identify  ferrous  objects", 
Proc.  IMACS  International  Symposium  on  A  I,  Expert  Systems  and  Languages  in  Modelling  and 
Simulation,  Institute  for  Cybernetics,  Barcelona,  Spain,  p.203,  Jun.4-6,  1987. 

[24]  D.Cohen,  “Developments  in  ways  to  measure  the  extremely  weak  magnetic  fields  emanating  from 
organs...",  Physics  Today,  pp.35-43,  August  1975. 

[25]  S.Breiner,  Applications  Manual  for  Portable  Magnetometers,  Sunnyvale  California:  GcoMetrics 
Ltd.,  1973. 


Inverse  Methods 

[26]  T.Caffey  and  L.Romero,  “Locating  a  buried  dipole",  IEEE  Trans  Geoscience  and  Remote  Sensing, 
GE-80,  pp.  188- 192,  1982. 

[27]  A.McAulay,  “Computerised  model  demonstrating  magnetic  submarine  localisation", IEEE  Trans 
AES,  AES-13,  pp.246-254,  1977. 

[28]  F.Raab,  E. Blood, T.Steiner  and  H.Jones,  “Magnetic  position  and  orientation  tracking  system",  IEEE 
Tram  AES,  AES-15,  pp.709-717,  1979. 

[29]  W.Wynn,  C.Frahm,  P.Carroll,  R.Clark,  J.  Weilhoner  and  M.  Wynn,  “Advanced  superconduct¬ 
ing  gradiometer/magnetometer  arrays  and  a  novel  signal  processing  technique",  IEEE  Trans  Mag, 
MAG-11,  pp.  701-707, 1975. 

[30]  J.McFee  and  Y.Das,  “Determination  of  the  parameters  of  a  dipole  by  measurement  of  its  magnetic 
field", IEEE  Trans  AP,  AP-B9 ,  pp.282-287,  1981. 

[31]  J.McFee  and  Y.Das,  “Fast  nonrecursive  method  for  estimating  location  and  dipole  moment  compo¬ 
nents  of  a  magnetic  dipole",  IEEE  Trans  GeoScience  and  Remote  Sensing,  GB-B4,  p.663,  1986. 

[32]  J.E.McFee  and  Y.Das,  “A  classifier  for  feature  vectors  whose  prototypes  are  a  function  of  multiple 
continuous  parameters",  IEEE  Trans.  PAMI,  1988. 

Measurement  Facilities 

[33]  J.McFee,  M.Bell,  B. Dempsey,  R.Chesney  and  Y.Das,  “A  magnetostatic  signature  measurement  and 
analysis  system" ,  J.Phys.E.,  18,  pp. 54-60,  1985. 


57 


Magnetometer  Types 


[34]  A.L.Bloom,  “Optical  pumping”,  Scientific  American ,  pp.72-80,  1860. 

[35]  P.Grivet  and  L.Malnar,  “Measurement  of  weak  magnetic  fields  by  magnetic  resonance",  Advances 
in  Electronics  and  Electron  Physics,  23,  pp.39-151, 1967. 

[36]  F.Primdahl,  “The  fluxgate  magnetometer”,  J.Phys.E.:Sci.Instrvm,  IS,  pp.241-253, 1979. 

[37]  F.Acker,  “Calculation  of  the  signal  voltage  induced  in  a  toroidal  proton  precession  magnetometer 
sensor”,  IEEE  Trans  Geoscience  Electronics,  GE-9,  pp.98-103,  1971. 

[38]  T.Casselman,  “Calculation  of  the  performance  of  a  magnetorestrictive  permalloy  magnetic  field 
sensor”,  IEEE  Trans  Magnetics,  MAG-16,  pp.461-464,  1980. 

[39]  J. Clarke,  “Low-  frequency  applications  of  superconducting  quantum  interference  devices”, 
Proc.IEEE  61,  pp.8-19, 1973.  stable  r.f.  biased  superconducting  point  -  contact  quantum  devices”, 
J.Appl.Phys.  41,  p.1572, 1970. 

[40]  A.Yariv  and  H.Winsor,  “Proposal  for  detection  of  magnetic  fields  through  magnetorestrictive  per¬ 
turbation  of  optical  fibers",  Optics  Letters,  5,  pp.87-89, 1980. 


58 


ELECTROSTATIC  METHODS 


General 

[41]  L.Slichter  and  M.Telkes,  “Electrical  properties  of  rocks  and  minerals”,  Handbook  of  Physical  Con¬ 
stants,  Geological  Society  of  America  Special  Paper,  F.Birch,  Ed.,  pp.299-319,  1942. 

[42]  E.Parkhomenko,  Electrical  Properties  of  Rocks .  New  York:  Plenum,  1967. 

[43]  M.J  .Campbell  and  J.Ulrichs,  “Electrical  properties  of  rocks  and  their  significance  for  lunar  radar 
observations”,  J.Geophys.Res.  74,  pp.  6867-  5881, 1969. 

[44]  D.Jupp  and  K.Vosoff,  “Stable  Iterative  Methoda  for  the  Inversion  of  Geophysical  Data",  Geophys. 
J.  R.  Asir.Soc.  42,  pp.957  -  976,  1975. 

Conductivity  Methods 

[45]  C.V.Keller  and  F.C.Frischknecht,  Electrical  Methods  in  Geophysical  Prospecting.  Oxford:  Pergamon, 
1966. 

[46]  G.Kunetz,  Principles  of  Direct  Current  Resistivity  Prospecting.  Berlin:  Gebruder  Borntrager,  1966. 


Conductivity  Imaging 

[47]  R.Henderson  and  J. Webster,  “An  impedance  camera  for  spatially  specific  measurements  of  the 
thorax",  IEEE  Trans. Bio. Eng.  BME-2S,  pp.260-264, 1978. 

[48]  R.J. Lytle  and  K.A.Dines,  “An  impedance  camera:  A  system  for  determining  the  spatial  varia¬ 
tion  of  electrical  conductivity”,  Report  UCRL-52413,  Lawrence  Livermore  Laboratory,  Livermore, 
California,  1978. 

[49]  K.A.Dines  and  R.J. Lytle,  “Computerised  geophysical  tomography",  Proc.IEEE  67,  p.1065,  1979. 

[50]  K.A.Dines  and  R.J.Lytle,  “Analysis  of  electrical  conductivity  imaging”,  Geophysics  43,  pp.  1025- 
1038, 1981. 

[51]  L.R.Price,  “Electrical  impedance  computed  tomography  (ICT):  A  new  imaging  technique”,  IEEE 
Trans.Nuc.Sci.  NS-26,  pp.2738-2739,  1979. 

[52]  R.Bates,  G. McKinnon  and  A.Seager,  “A  limitation  on  systems  for  imaging  electrical  conductivity 
distributions",  IEEE  Trans.Bio.Eng.  BME-2 7,  pp. 418-420,  1980. 

[53]  R.Nakonechny,  “A  Preconditioned  Conjugate  Gradient  Method  Using  a  Sparse  Linked-List  Tech¬ 
nique  for  the  Solution  of  Field  Problems",  M.Sc.  Thesis,  Technical  Report  TR8S-3,  University  of 
Manitoba,  1983. 

[54]  A.Wexler.B.Fry  and  M. Neuman,  “An  impedance  computed  tomography  algorithm  and  system", 
Proc.  Optical  Society  of  America  meeting  on  Industrial  Applications  of  Computed  Tomography  and 
NMR  Imaging,  Heda  Island,  Manitoba,  Canada,  August  1984. 

[55]  A.Wexler  and  C.J.Mandel,  “An  impedance  computed  tomography  algorithm  and  system  for  ground 
water  and  hasardous  waste  imaging”,  Proceedings  of  the  Second  Canadian/American  Conference  on 
Hydrogeology,  Banff  Alberta  June  25-20  1985.  (pub.  National  Water  Well  Assoc.,  Dublin  Oil) 


59 


[56]  A.Tamburi,R.  Allard  and  U.Roeper,  “Tomographic  imaging  of  ground  water  pollution  plumes”,  Pro¬ 
ceedings  of  the  Second  Canadian/ American  Conference  on  Hydrogeology,  Banff  Alberta  June  25-29 
1985.  (pub.  National  Water  Well  Assoc.,  Dublin  OH) 


Induced  Polarisation 

[57]  J.S. Sumner,  Principles  of  Induced  Polarization  for  Geophysical  Exploration.  Amsterdam:  Elsevier, 
1976. 

[58]  -,  Advances  in  Induced  Polarization  and  Complex  Resistivity  (short  course ).  Tucson:  University  of 
Arisona,  1981. 

[59]  A.Howland-Rose,  G. Linford,  D. Pitcher  and  H.Seigel,  “Some  recent  magnetic  polarisation  develop¬ 
ments  -  Part  l:theory”,  Geophysics  45,  pp.37-  54,  1980. 


60 


ELECTROMAGNETIC  INDUCTION 


General 

[60]  D.F.Bleil,  Ed.,  Natural  Electromagnetic  Phenomena  Below  30  KC/S ,  New  York:  Plenum,  1964. 

[61]  J  E.Hipp,  “Soil  electromagnetic  parameters  as  functions  of  frequency,  soil  density  and  soil  moisture”, 
*  '.IEEE  62,  p.88,  1974. 


Electromagnetic  Induction  Theory 

[62]  J.R.Wait  and  K.P.Spies,  “Quasi-static  transient  response  of  a  conducting  permeable  sphere”,  Geo¬ 
physics  35,  pp.789-  792,  1969. 

[63]  M.N.Nabighian,  “  Quasi-static  transient  response  of  a  conducting  permeable  sphere  in  a  dipolar 
field”,  Geophysics  35,  pp.303-  309,  1970. 

[64]  S.K.Singh,  “Electromagnetic  transient  response  of  a  conducting  sphere  embedded  in  a  conducting 
medium",  Geophysics  38,  pp.864-  893,  1973. 

[65]  J.R.Wait,  “A  conducting  permeable  sphere  in  the  presence  of  a  coil  carrying  an  oscillating  current”, 
Con.  J.  Phys.  31,  pp.  670-  678,  1953. 

[66]  Y.Das,J.E.McFee  and  R.H.Chesney,  “Time  domain  response  of  a  sphere  in  the  field  of  a  coil-theory 
and  experiment",  IEEE  Trans  GeoScience  and  Remote  Sensing  GE-22,pMQ,  1984. 

[67]  K.M.Mitsner,  “An  integral  equation  approach  to  scattering  from  a  body  of  finite  conductivity", 
Rad.Sci.  V2,  1967. 

[68]  T.B.A.  Senior,  “Impedance  boundary  condition  for  imperfectly  conducting  surfaces",  App.Sci.Res. 
8B,  1961. 

[69]  A.A.Sebak  and  L.Shafai,  “Transient  response  computation  of  spheroidal  objects  using  impedance 
boundary  conditions", IEEE  Trans  AP  AP-32,  p.1118,  1984. 

Applications 

[70]  J.R.Wait,  “Review  of  electromagnetic  methods  in  nondestructive  testing  of  wire  ropes",  Proc.IEEE. 
67, pp.892-  903,  1970. 

[71]  R.Callarotti  and  M.  Alfonso,  “Measurement  of  the  conductivity  of  metallic  cylinders  by  means  of  an 
inductive  method",  J.  Appl.  ,Phys.  43,  pp.  3040-3047,  1972. 

[72]  C.Colani  and  M.Aitken,  “A  new  type  of  locating  device.  The  instrument",  Pwsspesioni  Archco- 
logiche,  pp.15-23,  I960. 

[73]  E.B.Glennie,  “Inductive  detection  of  buried  services",  Proe.IEB  125,  pp.14-15,  1978. 


Inverse  Methods 

[74]  Y.Das,  J.E.McFee  and  R.H.Chesney,  “Determination  of  depth  of  shallowly  buried  objects  by  elec¬ 
tromagnetic  induction", ,/.£££  Trans  GeoScience  and  Remote  Sensing  GE-2S,  p.00,  1985. 

[76]  S.H.Ward,  “Unique  determination  of  conductivity,  susceptibility,  sise  and  depth  in  multifrequency 
electromagnetic  exploration",  Geophysics  U,  p.531,  1959. 


61 


[76]  W.H.Huggins,  "Representation  and  analysis  of  signals  part  1,  the  use  of  orthogonalised  exponen¬ 
tials”,  Department  of  Electrical  Engineering  Report,  Johns  Hopkins  University,  1957. 

[77]  R.N.McDonough,  "Matched  exponents  for  the  representation  of  signals”,  Ph.D. Thesis  Johns  Hop¬ 
kins  University,  1963. 

[78]  Y.Daa,J.E.McFee  and  M.Bell,  “Detection  and  identification  of  buried  ordance  by  magnetic  and 
electromagnetic  means  (U)”,  Suffield  Report  283  1981,  UNCLASSIFIED. 

[79]  R.H.Chesney,  Y.Das,  J.E.McFee  and  M.Ito,  "Identification  of  metallic  spheroids  by  classification  of 
their  electromagnetic  responses”,  IEEE  Trans  PAMI PAMI-6, p.809,  1984. 


Instruments  and  Facilities 

[80]  B.Roston,  “Development  of  locators  of  small  metallic  bodies  buried  in  the  ground”,  J.  of  the  IEE 
95,  pp.653-  667,  1948. 

[81]  R.Macario,  “Discriminative  metal  detector”,  Wireless  World,  pp.43-48,  July  1978. 

[82]  J.Corbyn,  “Pulse  induction  metal  detector”,  Wireless  World,  pp.40-44,  March  1980. 

[83]  J.Corbyn,  “Pulse  induction  metal  detector  -  2”,  Wireless  World,  pp.58-59,  April  1980. 

[84]  J.E.McFee,  R.H.Chesney,  Y.Das,  J.Toews,  M.Turnbull,  R.Pinnell,  and  M.Bell,  “An  experimental 
time  domain  electromagnetic  induction  system”, Rev.  Sci.  Instrum  55,  p.968,  1984. 

[85]  Y.Das,  J.Toews  and  J.E.McFee,  “Vehicle-mounted  ordnance  locator  :  an  experimental  prototype 
(U)”,  Suffield  Report  431  1988,  UNCLASSIFIED. 


62 


MAIN  CHARACTERISTICS  OF  EXPLOSIVE  MUNITIONS 
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discrimination 
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Table  1.2 

METHODS  FOR  DETECTING  EXPLOSIVE  MUNITIONS 


METHOD 


REMOTE/ 
QUASI -REMOTE 


Magnetostatics  Q 

Electromagnetic  Induction  Q 

Conductivity  Imaging  Q 

Ground  Penetrating  Radar  Q 

Acoustics  Q 

RF  Resonance  Absorption  Q 

Nuclear  Q 

Trace  Gas  Analysis  Q 

Biochemical  Detection  Q 

Optical  -  Ultraviolet  Q,,  R 

Optical  -  Visible  R 

Optical  -  Infrared  -  Near  R 

Optical  -  Infrared  -  Thermal  Q,  R 

Microwaves  ~  Passive  Q,  R 

Microwaves  -  Active  Q,  R 

Human  Perception  Channels  ? 


UNCLASSIFIED 
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Table  1.3 

ELECTRICAL  PROPERTIES  OF  TWO  COMMON 
MATERIALS  IN  REMOTE  SENSING 


MATERIAL 

CONDUCTIVITY 

FREQUENCY 

SKIN  DEPTH 

S/m 

Hz 

m 

SALINE 

101 

1000 

0.16 

100 

0.60 

10 

0.16 

SOIL 

10“ 1 

1000 

1.6 

(MAXIMUM 

100 

5.0 

CONDUCTIVITY) 

10 

16. 

UNCLASSIFIED 
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Table  2.1 

STANDARD  NOMENCLATURE  OF 
GEOMAGNETIC  FIELD  COMPONENTS 


Component 

Standard 

letter 

symbol 

Average  value  near 
latitude  45°N 
(nT) 

Total  field 

F 

45000 

Horizontal  component 

H 

15000-25000 

Vertical  component 

Z 

40000 

Declination  (angle  between  H  and  the 
geographical  north  direction) 

D 

Depends  on  location 

Inclination  or  dip  (angle  between  F 
and  the  horizontal  plane) 

I 

60-70°  (F  pointing  down) 

UNCLASSIFIED 
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Table  2.2 

MAGNETIC  UNITS 


Quantity 

Typical 

Symbol 

Units 

Magnetic  Induction  or 

6,  B 

1  nanoTesla  (nT) 

Magnetic  Flux  Density 

=  10~9  Tesla  (Weber/m2) 

—  10~5  Gauss 
=  1  gamma  (-7) 

Magnetic  Field  or 

hjj 

1  uT/mo 

Magnetic  Field  Intensity 

=  10'5  Oersteds  (Gilberts/cm) 

=  7.00  x  10""  A/m 
=  7.96  x  10”7  emu 

Magnetization 

M' 

same  as  H 

Magnetic  Dipole 

m,  M,  Aft1) 

1  Am3 

Moment 

-  103emu 

UNCLASSIFIED 


COMPARISON  OF  SOME  PROPERTIES  OF  Q  AND  B3  IN  A 
HORIZONTAL  MEASUREMENT  PLANE  (See  Figure  2.13) 
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UNCLASSIFIED 


Determination,  of  depth  not  straightforward  comes  from  same  computation 

and  dipole  moment  from  simple  measurements  as  for  position  in  plane 
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Table  2.4 

DIMENSIONS  OF  SPHEROIDAL  OBJECTS  FOR 
CLASSIFICATION  EXPERIMENTS 


OBJECT 

NUMBER 

LENGTH  OF  SYMMETRY 
AXIS  (a)  IN  M 

RATIO  OF  AXES  LENGTHS 

M _ 

1 

0.010 

3.75 

2 

0.020 

3.40 

3 

0.030 

2.50 

4 

0.045 

2.60 

5 

0.060 

3.50 

8 

0.090 

3.50 

UNCLASSIFIED 
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Table  2.5 

PROBABILITY  OF  MISCLASSIFICATION  FOR 
SPHEROIDS  OF  (  ABLE  24 


PERCENT 

NOISE 

PROBABILITY  OF 
MISCLASSIFICATION  (%) 

PERCENT 

NOISE 

PROBABILITY  OF 
MISCLASSIFICATION  (%) 

0 

1.1 

10 

7.1 

1 

1.1 

25 

16.4 

2 

1.4 

50 

30.3 

5 

2.9 

75 

39.0 

UNCLASSIFIED 


UNCLASSIFIED 
Table  3.1 
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ELECTROMAGNETIC  PROPERTIES  OF 
RELEVANT  MATERIALS 


M&terial 

Dielectric 

Conductivity 

Relative 

Constant 

(S/m) 

Permeability 

A.Soils 


dry  sand 

4  — ►  6 

10-7  -*•  10”3 

1 

water  saturated  sand 

30 

10”4  -¥  IQ~2 

1 

water  saturated  silt 

10 

10”3  10”3 

1 

water  saturated  clay 

8  -*  12 

10”1  ->  1 

1 

dry, sandy, flat  coastal  land 

10 

2  X  10”3 

1 

marshy, forested, flat  land 

12 

8  x  10-3 

1 

rich  agricultural  land,  low  hills 

15 

10”  3 

1 

pastoral  land,  medium 
hills  and  forestation 

13 

5  x  10”3 

1 

granite,  dry 

5 

10”® 

1 

limestone,  dry 

7 

10”® 

1 

B.Metate 


silver 

1 

6,3  x  10+7 

1 

1 

6,0  X  10+7 

1 

gold 

1 

4.3  X  10+7 

1 

aluminum 

1 

3.8  x  10+7 

1 

iron 

1 

1.0  x  10+T 

50  —*  1000 

bismuth 

1 

9.4  x  10+6 

1 

C. Miscellaneous 

air 

1 

0 

fresh  water 

81 

10”4  -*3x  10”3  1 

sea  water 

81 

4 

fresh  water  icc 

4 

10”4  ->  10”3  1 

permafrost 

4  — ♦  8 

10”s  -*  10 *3  1 

UNCLASSIFIED 
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Table  3.2 

CONDUCTIVITIES  FOR  THREE  SIMPLE 
ANISOTROPY  MODELS 


Model 

o\  =  lOoj 

0\  —  lOOoj 

a 

at  =  oy  =  a,  =  1.47ctj 

ot  —  Oy  —  oM  —  7A7oj 

b 

oz  —  ov  —  2.8<7j 

Oz  —  Oy  —  20.8O] 

a,  =  1.22aj 

ot  =  1.250} 

c 

oz  —  oy  —  ot  —  X.53<Tj 

O*  =  Oy  -  O,  =  1.72(72 

UNCLASSIFIED 
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MAGNETIC 

NORTH 


GEOGRAPHIC 

NORTH 


Rgure  2.1 

THE  EARTH  AS  A  MAGNETIC  DIPOLE 
(Grivet  and  Malnar) 


UNCLASSIFIED 
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UNCLASSIFIED 


Figure  2.2 

THE  GEOMAGNETIC  INCLINATION  IN  DEGREES  OF  ARC  FROM  THE  HORIZONTAL 


135°  180°  135° 
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Figure  2.3 

THE  TOTAL  INTENSITY  OF  THE  EARTH'S  MAGNETIC  FIELD 
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MID-NORTHERN  AND  MID -SOUTHERN  LATITUDES 


TYPICAL  DIURNAL  VARIATIONS  IN  TOTAL  FIELD  INTENSITY 

Figure  2.4 


Figure  2,8 


TYPICAL  MAGNETIC  STORM 
Figure  2.6 

TEMPORAL  VARIATIONS  IN  THE  GEOMAGNETIC  FIELD  (Breiner) 

UNCLASSIFIED 
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Figure  2.7 

PIONEER  V  "NEAR  ZONE'  FIELD  MEASUREMENTS 
l  G rivet  and  Malnarj 
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MAGNETIC  INDUCTION  {NT)  TYPICAL  SOURCES 
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Figure  2.8 

TYPICAL  MAGNETIC  SOURCE  STRENGTHS  AND  DETECTOR  SENSITIVITIES 
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Figure  2.9 

FLUXGATE  MAGNETOMETER  PRINCIPLE  OF  OPERATION 
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EARTH'S 

MAGNETIC  FIELD 
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SENSORS  ALIGNED 


FIELD  FROM : 


60000 

60000 

EARTH 

59088 

59088 

200 

100 

+  OBJECT 

197 

98 

60200 

60100 

TOTAL 

59285 

59186 

MEASURED 


a  2  *  F,ELD  b 


1 

2 

i 

2 

-B2  +  B, 

-  100  nT 

-b2  +  B, 

-  99  nT 

SENSORS  MISALIGNED  1° 

60000 

200 

59991 

100 

59088 

197 

58897 

98 

60200 

60091 

59285 

58995 

2 

2 

d 

-B2  +  B,  a  109  nT  -B2  +  B,  =  290  nT 

Figure  2.10 

HOW  MISALIGNMENT  IN  A  VECTOR  GRADIOMETER  AFFECTS 
ERROR  IN  GRADIENT  MEASUREMENT 

UNCLASSIFIED 


CART  WITH 
MAGNETOMETER 
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POSITION  ON  1  AXIS  OF  MAXIMUM  OF  THE  GRADIENT  MEASUREMENT  QUANTITY  Q 
(NORMALIZED  TO  DEPTH  OF  DIPOLE)  VERSUS  RATIO  OF  1-  AND  3-  COMPONENTS  OF 

DIPOLE  MOMENT 


UNCLASSIFIED 
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nT 


4w 


W,er*BOVE  A  0*0^ 

— ,A‘m 
Dipole  Moment  <Wv  "2 


unclassified 


figure  2.1® 

,mv  ‘ 


x  0»POl£ 
.  0^  &,rri 


UNCtASSlP'^0 


ft ne  i ", 
rn2.  m3>  “ 


0.  a)  A  t" 


^classified 


VALUE  OF  1/2  Qmax  *«  PLANE  1  m  ABOVE  A  DIPOLE 
Note  that  focus  is  not  centered  about  ( X-j,  X2)  =  (0,  0). 


GEOMETRY  FOR  CALCULATION  OF  MULTIPOLE  MOMENTS  OF  A  SPHEROID 
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SPHEROID  SHAPE  PARAMETERS  e 
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RATIO  OF  Fj  TO  F f  VERSUS  SPHEROID  RELATIVE  PERMEABILITY  FOR  VARIOUS 

SPHEROID  SHAPE  PARAMETERS  a 

Relative  permeability  of  surrounding  medium  Is  1.  For  fixed  orientation,  ratio  will  be 
proportional  to  ratio  of  dipole  moments  induced  in  spheroid  along  3  and  1  axes. 


f3/f, 


UNCLASSIFIED 


Figure  2.26 

RATIO  OF  F.  TO  F,  VERSUS  RELATIVE  PERMEABILITY  OF  SURROUNDING  MEDIUM  FOR 

VARIOUS  SPHEROID  SHAPE  PARAMETERS  o 
Relative  permeability  Of  spheroid  Is  SO.  Most  rocks  have  pr2  -  1  vvhile  the 
maximum  for  most  ferrous  rocks  is  (if2  ~  3’ 


UNCLASSIFIED 
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POSITION  FROM  ORIGIN  (m) 

Figure  2.27 

MEASURED  TOTAL  MAGNETIC  FIELD  OF  A  .08  m  RADIUS  MILD  STEEL  SPHERE 

AT  A  DEPTH  OF  0.73m  (SOLID  LINE) 

Geometry  is  that  of  Figure  2.11.  Sphere  center  is  directly  under  position 
value  0.  Dotted  line  is  theoretical  prediction  (Me Fee  et  aJ„  1985). 


c 


POSITION  FROM  ORIGIN  !m) 
Figure  2.28 


MEASURED  TOTAL  MAGNETIC  FIELD  OF  A  105  mm  HJWiT ZER  SHELL  AT  A 

DEPTH  OF  0.96m  (SOLID  LINE). 

Geometry  is  that  of  Figure  2.11.  Symmetry  axis  is  tilted  47°  from  vertical  and 
projection  is  43°  from  magnetic  north.  Shefl  center  is  directly  under  position  value 
Dotted  line  is  theoretical  prediction  (Me Fee  et  a/.,  1985). 


UNCLASSIFIED 


P  ff 


UNCLASSIFIED 


SSP  124 


0  =  0 


Figure  2.29 

THE  DESIGN  PROTOTYPE  MANIFOLDS  FOR  TWO  DIFFERENT  FERROUS  SPHEROIDS. 
Each  manifold  is  a  surface  consisting  of  a  continuum  of  points.  Each  point  is  the 
head  of  a  magnetic  dipole  moment  vector  corresponding  to  a  polar  angle,  0,  and 
a2imuthai  angle,  of  the  symmetry  axis  of  the  sphoroid  with  respect  to  a 

space-fixed  coordinate  system. 

UNCLASSIFIED 


UNCLASSIFIED 


SSP  124 


Figure  2.30 

PORTION  OF  A  TWO  DIMENSIONAL  PROTOTYPE  MANIFOLD  IN  A 
THREE  DIMENSIONAL  FEATURE  SPACE  AND  ITS  APPROXIMATION 

BY  A  SET  OF  UNIT  CELLS. 

Each  unit  cell  Is  a  triangle  formed  by  three  adjacent  points  sampled 

on  the  manifold. 
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CO ND U CTtVJTY  SPECTRUM  (Ward  and  Fraser) 
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Figure  32 

SCHEMATIC  REPRESENTATION  OF  THE  WENNER  ARRAY 
The  parameter  a  would  be  increased  to  measure  deeper  layers  in  the  medium. 
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Figure  3.3 

TWO  SCHEMES  FOR  HORIZONTAL  EXPLORATION 
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Figure  3.4 

THREE  MODELS  OF  CONDUCTIVITY  ANISOTROPY 
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IGNEOUS 

ROCKS 


METAMORPHIC 

ROCKS 


SEDIMENTARY 

ROCKS 


UNCONSOLIDATED 
SEDIMENTS 
AND  SOILS 


CONDUCTIVITY, 

S/m 

RESISTIVITY, 


ohm-rn 


Figure  3.5 

CONDUCTIVITIES  OF  TYPICAL  ROCK  TYPES 
(Slichter  and  Telkes) 


UNCLASSIFIED 


UNCLASSIFIED 


SSP  124 


UNCLASSIFIED 
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Figure  3.7 

TWO  LAYER  EARTH  GEOMETRY 
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Figure  3.8 

BLOWUP  OF  THE  ELECTRODE  PORTION  (0,0,0)  OF  FIGURE  3.7 
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Figure  3.9 

FOUR  ELECTRODE  GEOMETRY 


UNCLASSIFIED 


UNCLASSIFIED 


SSP  124 


Figure  3.10 

RATIO  OF  APPARENT  RESISTIVITY  TO  TOP  LAYER  RESISTIVITY  FOR  2  LAYER 
EARTH  MODEL  AND  WENNER  ARRAY  OF  SPACING  a  (Wait  1987.) 

Top  layer  thickness  is  h .  K  is  defined  In  text 
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CURRENT 


Figure  3.11 

GEOMETRY  FOR  CONDUCTIVITY  IMAGING  OF  HIDDEN  OR  BURIED  COMPACT  OBJECTS 
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UNCLASSIFIED 


DISTRIBUTION  BASED  ON  SIMULATED  MEASUREMENT  OF  FIELDS  ON  THE 

TOP  SURFACE  OF  THE  VOLUME  (Wexler) 
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Figure  3.13 

MICROSCOPIC  MODEL  FOR  INDUCED  POLARIZATION  IMPEDANCE 
Spherical  particle  has  a  complex  surface  impedance  n  and  volume  resistivity  q 
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Figure  3.14 

MACROSCOPIC  SHPERICAL  VOLUME  CONTAINING  A  NUMBER  OF  IDENTICAL 
MICROSCOPIC  SPHERES  AS  IN  FIGURE  3.13 
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Figure  3.15 

NORMALIZED  STEP  FREQUENCY  RESPONSE  USING  THE  W-POLE 

MODEL  FOR  RESISTIVITY 


Figure  3.16 


TYPICAL  IMPRESSED  CURRENT  WAVEFORM  USED  IN  INDUCED 
POLARIZATION  EXPERIMENTS 


Figure  3.17 


UNNORMALIZED  RESPONSE  TO  WAVEFORM  OF  FIGURE  3.16  USING 
THE  W-POLE  MODEL  FOR  RESISTIVITY 
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FREQUENCY  (Hz) 


THE  GEOMAGNETIC  FREQUENCY  SPECTRUM  (Bteil) 
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Figure  4.4 

GEOMETRY  FOR  MUTUAL  INDUCTANCE  CALCULATIONS  FOR 
TWO  CURRENT  CARRYING  LOOPS 
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Figure  4.5a.  COAXIAL  GEOMETRY 
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Figure  4.3b.  COPLANAR  GEOMETRY 


Figure  43c.  ORTHuGONAL  GEOMETRY 


Figure  4.6 

SPECIAL  CASES  FOR  MUTUAL  INDUCTANCE  CALCULATIONS 

For  all  cases  f  t>avj2 
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Figaro  4.6a 
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Figure  4.6  c 
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Figure  4.6 

VARIOUS  COIL  COUPLING  GEOMETRIES  FOR  USE 
OVER  AN  INFINITE  FLAT  SHEET 
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Figure  4.7 

SIMPLE  CIRCUIT  MODEL  OF  ELECTROMAGNETIC  INDUCTION 
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INDUCTIVE  LIMIT  (REAL) 
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INDUCTION  (Grant  and  West) 
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PHASOR  DIAGRAMS  SHOWING  EVOLUTION  OF  INDUCED  CURRENT  IN 

TARGET  CIRCUIT  AS  FREQUENCY  G)  OF  INDUCING  FIELD  INCREASES 

* 

V1  is  due  to  primary  magnetic  field  (current  Iq).  V1  is  due  to  secondary 

field  (current  ly )  ( Grant  and  West). 
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TYPICAL  RESPONSE 
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RESPONSE  AS  A  FUNCTION  OF  POSITION  OF  SENSOR  SYSTEM  WITH 

RESPECT  TO  TARGET  POSITION 
Geometry  is  that  of  Figures  4.7,  4.11  and  4.12. 
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TRANSMITTER  RECEIVER  TRANSMITTER  RECEIVER 


Figure  4.11 

FLUX  DIAGRAM  TO  QUALITATIVELY  EXPLAIN  RESPONSE  OF  FIGURE  4.10 
In  this  case  transmitter/receiver  straddle  target. 


TRANSMITTER  RECEIVER  TRANSMITTER  RECEIVER 


Figure  4.12 

FLUX  DIAGRAM  TO  QUALITATIVELY  EXPLAIN  RESPONSE  OF  FIGURE  4.10 
In  this  case  transmitter /receiver  are  to  one  side  of  target 
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Z  PRIMARY  LOOP 


Figure  4.13 

GEOMETRY  FOR  CALCULATION  OF  THE  RESPONSE  OF  A 
SPHERE  IN  THE  FIELD  OF  A  COIL 
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RESPONSE  FUNCTION  (X„  +  /yn) 
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RESPONSE  PARAMETER  (a/i0O)a2) 


Figure  4.14 

RESPONSE  FUNCTIONS  OF  INDUCED  MULTIPOLE  MOMENTS  FOR  A 
SPHERICAL  CONDUCTOR  IN  THE  FIELD  OF  A  COIL 
The  4  lowest  orders  are  shown  (Wait  1953). 
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FIRST  10  TERMS  OF  THE  TIME  RESPONSE  OF  A  MAGNETICALLY 
NONPERMEABLE  SPHERE  IN  A  UNIFORM  MAGNETIC  FIELD 
(Das  and  McFoe,  1981 K 


FIRST  10  TERMS  OF  THE  TIME  RESPONSE  OF  A  MAGNETICALLY 
PERMEABLE  SPHERE  IN  A  UNIFORM  MAGNETIC  FIELD 
( Das  and  McFse,  1981 ). 
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FEATURE  SPACE  COORDINATE  1 
Figure  4.18 

PORTION  OF  A  ONE  DIMENSIONAL  PROTOTYPE  MANIFOLD  (A  CURVE) 
IN  AN  N  DIMENSIONAL  FEATURE  SPACE 
(To  allow  it  to  be  illustrated,  the  manifold  has  been  projected  onto  the  2-space 
formed  by  two  of  Its  feature  coordinates).  Its  epproximetlon  by  a  set  of 

line  segment  unit  cells  is  shown. 

UNCLASSIFIED 
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